ERP model
noslopes=TRUE # set to use just random intercept model
noslopes=FALSE # set to false to get maximal model
if (!file.exists("resultdata.rds")){
write("reading dataframe.csv", stderr())
# load data for all models
# resultWhole.df <- read.csv("smalldevel.csv")
resultWhole.df <- read.csv("dataframe.csv") # automatically uses working dir
write("finished reading dataframe.csv", stderr())
result.df <- resultWhole.df[,0:97]
# result.df <- removeNAColumns(resultWhole.df)
names(result.df)[6]<-"Condition"
result.df$file=as.character(result.df$file)
result.df$sub = str_extract(result.df$file,"[-]s[0-9]+[-]")
result.df$sub = str_replace_all(result.df$sub,"[-]","")
result.df$sub=factor(result.df$sub)
result.df$epoch = as.integer(str_extract(result.df$file,"[0-9]+$"))
result.df$coul = str_extract(result.df$file,"[0-9.]+vcoul[12]")
result.df$vrom = str_extract(result.df$file,"[0-9.]+vrommer")
result.df$vromNot = str_extract(result.df$file,"[0-9.]+vromNot")
result.df$devellearn = str_extract(result.df$file,"[0-9.]+violation")
# result.df$file = NULL
epochCol = which(names(result.df)=="epoch")
subCol = which(names(result.df)=="sub")
coulCol = which(names(result.df)=="coul")
vromCol = which(names(result.df)=="vrom")
vromNotCol = which(names(result.df)=="vromNot")
devellearnCol = which(names(result.df)=="devellearn")
result2.df = result.df[,c(subCol,coulCol,vromCol,vromNotCol,devellearnCol,epochCol,0:(subCol-1))]
names(result2.df)[which(names(result2.df)=="sent")]<-"Example"
result2.df$tick = result2.df$tick+1
result2.df$layer = factor(result2.df$layer, levels=c("word","compress", "hidden","cword"))
result2.df$depth = as.numeric(result2.df$layer)
result.df=result2.df
# saveRDS(result.df,file="smalldevel.rds")
saveRDS(result.df,file="resultdata.rds")
}else{
print("using Rds file")
# devresult.df=readRDS(file="smalldevel.rds")
result.df=readRDS(file="resultdata.rds")
}## [1] "using Rds file"#print folder list to make sure correct folders
folder = str_replace(result.df$file[1],"/derivatives.*","")
if (folder == "derivatives100000"){
folder = basename(getwd())
}
sublist = paste(levels(result.df$sub),collapse=",")
print(paste("Folder",folder," nsub=",sublist))## [1] "Folder sim80-s0-Sim116-alt10-50/deriv0.1vcoul1-100000 nsub= s0,s1,s2,s3,s4,s5,s6,s7,s8,s9"# load counts
counts.df = read.csv("counts.csv")
counts.df$perc = round(100*counts.df$prop)
#print(counts.df[order(counts.df$pair),])
unigram.df = read.csv("unigram.csv")
#print(unigram.df[unigram.df$pair%in%c("wine","water","tea","cake"),])
result2.df=result.df
first=which(names(result2.df)=="the")-1
end=which(names(result2.df)=="per")
# compute metrics Hidden layer is smaller than Lexical layer, so this captures all units
# hidden has only 50 units, so rest are marked with NA
result2.df$abssum = rowSums(abs(result2.df[,first:end]),na.rm = TRUE) # sum abs error
result2.df$absmean = rowMeans(abs(result2.df[,first:end]),na.rm = TRUE) # mean abs error
object.size(result2.df[,first:end])## 3793881248 bytesSys.getenv("R_MAX_VSIZE")## [1] "40Gb"couldf = subset(result2.df,epoch == max(result2.df$epoch) & !layer %in% c("cword") & !is.na(coul))
couldf$file = NULL
couldf$vrom = NULL
couldf$vromNot = NULL
#xtabs(~ Condition + coul, couldf)
print(unique(couldf$coul))## [1] "0.1vcoul1" "0.1vcoul2"rommerdf = subset(result2.df,epoch == max(result2.df$epoch) & !layer %in% c("cword") & !is.na(vrom))
rommerdf$file = NULL
rommerdf$coul = NULL
rommerdf$vromNot = NULL
rommerNotdf = subset(result2.df,epoch == max(result2.df$epoch) & !layer %in% c("cword") & !is.na(vromNot))
rommerNotdf$file = NULL
rommerNotdf$coul = NULL
rommerNotdf$vrom = NULL
develop.df = subset(result2.df,!layer %in% c("cword") & is.na(coul) & is.na(vrom) & is.na(vromNot))
develop.df$file = NULL
develop.df$coul = NULL
develop.df$vrom = NULL
develop.df$vromNot = NULL
develop.df$devellearn = NULL
adultdata.df = subset(develop.df,epoch == max(develop.df$epoch) )
adultdata.df$file = NULL
adultdata.df$vrom = NULL
adultdata.df$coul = NULL
adultdata.df$vromNot = NULL
adultdata.df$devellearn = NULL
#load("adultdata.RData")
layerNotERP = c("output","target", "error","myoutputder","myinputder")
dependMeasure = "abssum" # we can change the dependent measure for all models
mainMeasure = "inputderiv"
p600layer = "hidden"
#load(adultdata.df)
first=which(names(adultdata.df)=="the")
end=which(names(adultdata.df)=="per")
# when doing final version, change noslopes to FALSE
figwidth=6
figheight=3Cloze
write("cloze", stderr())
cloze.df = subset(adultdata.df, str_detect(adultdata.df$Condition,"CLOZE") & wordcat %in% c('NOUNI'))
cloze.df$Condition = factor(cloze.df$Condition,labels=c("High Cloze","Medium Cloze","Low Cloze"))
print(head(cloze.df[cloze.df$measure=="target",c("word","wordcat","Condition","Example","tick")])) # check## word wordcat Condition Example
## 114921 water NOUNI High Cloze a teacher was drink -ing the water . .
## 115038 water NOUNI Medium Cloze a teacher was taste -ing the water . .
## 115155 water NOUNI Low Cloze a teacher was take -ing the water . .
## 117053 water NOUNI High Cloze the driver will drink the water . .
## 117157 water NOUNI Medium Cloze the driver will taste the water . .
## 117261 water NOUNI Low Cloze the driver will take the water . .
## tick
## 114921 7
## 115038 7
## 115155 7
## 117053 6
## 117157 6
## 117261 6# check that layers are the right size
cloze.df$laylen = 1+(end-first)-rowSums(apply(cloze.df[,first:end],2,is.na))
print(aggregate(laylen ~ layer + measure, cloze.df,mean))## layer measure laylen
## 1 word inputderiv 87
## 2 compress inputderiv 29
## 3 hidden inputderiv 49
## 4 word outderiv 87
## 5 compress outderiv 29
## 6 hidden outderiv 49
## 7 word output 87
## 8 compress output 29
## 9 hidden output 49
## 10 word target 87cloze.df$laylen=NULL
showWordOutError("coffee,tea,wine,water,cake",c("Condition","measure"),cloze.df)
computeClozeERPCorr(cloze.df,dependMeasure,"word")## [1] "Correlation output prediction and abssum error of word is -0.999999999999812"computeClozeERPCorr(cloze.df,dependMeasure,"hidden")## [1] "Correlation output prediction and abssum error of hidden is -0.446316565145821"corrHiddenOutputInputDeriv(cloze.df,"Condition",layername="word")## Condition
## High Cloze Medium Cloze Low Cloze
## 8 8 8
myggsave("img/actderivcloze.png",width=6,height=6)# this is for creating ERP figure
nooutput = subset(cloze.df, ! measure %in% layerNotERP )
end = which(names(nooutput)=="depth")-1
first = which(names(nooutput)=="none")
norawcloze.df = nooutput[,-(first:end)] # remove word specific activations/error
dvlen = length(names(norawcloze.df))- which(names(norawcloze.df)=="depth")
longnoraw= melt(norawcloze.df,id.var=1:(length(norawcloze.df)-dvlen)) # make long format
print(folder) # we use folder rename html file## [1] "sim80-s0-Sim116-alt10-50/deriv0.1vcoul1-100000"# add target layer for figure
#mapping = aes(x=depth, y=value, colour = Condition,linetype=Condition)
p1 = drawERP(longnoraw,"value", "Condition" ,"measure","variable" ,timesize=-1)
#p1= drawERP(longnoraw,mapping, timesize=-1,span = 0.8)
p1 = p1 + facet_wrap( ~ measure+variable,scales = "free",nrow=1)
p1
print(folder)## [1] "sim80-s0-Sim116-alt10-50/deriv0.1vcoul1-100000"# this is the counts for verb arguments in the input
interactionText="NO TEXT"
# get example for table
tablesub = subset(norawcloze.df,layer == "word" & measure == mainMeasure)
exampleTable = tablesub[1:3,c("Condition","Example")]
tabletext = kable(exampleTable,format="latex",row.names=FALSE)
print(exampleTable)## Condition Example
## 114923 High Cloze a teacher was drink -ing the water . .
## 115040 Medium Cloze a teacher was taste -ing the water . .
## 115157 Low Cloze a teacher was take -ing the water . .wordccompinput = subset(norawcloze.df,measure == mainMeasure & layer %in% c("word",p600layer) )
if (nlevels(adultdata.df$sub) > 4){ # do mixed with more than 9 subj
print("run models")
# our analysis uses this set
wordccompinput$cloze = 0
wordccompinput$cloze[wordccompinput$Condition=="Low Cloze" ]=1
wordccompinput$cloze[wordccompinput$Condition=="High Cloze" ]=-1
wordccompinput$cword = ifelse(wordccompinput$layer=="word",0.5,-0.5)
# this is the main mixed model with centered variables
randomeff = " + (1 + cloze*cword |sub)"
if (noslopes){ randomeff = " + (1 | sub)" }
formu = as.formula(paste(dependMeasure,"~ cloze*cword",randomeff ))
omnimodel = lmer(formu, wordccompinput,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodel))
modellist = modelComparison(omnimodel) # create anova table by terms subtraction
}## [1] "run models"
## Linear mixed model fit by REML ['lmerMod']
## Formula: abssum ~ cloze * cword + (1 + cloze * cword | sub)
## Data: wordccompinput
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: 826.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8820 -0.4567 -0.1710 0.5325 6.0361
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.01018 0.1009
## cloze 0.03257 0.1805 -0.23
## cword 0.03299 0.1816 -0.47 -0.68
## cloze:cword 0.04520 0.2126 -0.30 -0.82 0.78
## Residual 0.08734 0.2955
## Number of obs: 1800, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.17013 0.03265 35.838
## cloze 0.41329 0.05771 7.162
## cword 0.67599 0.05910 11.438
## cloze:cword 0.19412 0.06936 2.799
##
## Correlation of Fixed Effects:
## (Intr) cloze cword
## cloze -0.222
## cword -0.442 -0.649
## cloze:cword -0.289 -0.786 0.733
## [1] "remove . ~ . - cloze:cword "## refitting model(s) with ML (instead of REML)## Data: wordccompinput
## Models:
## model2: abssum ~ cloze + cword + (1 + cloze * cword | sub)
## model: abssum ~ cloze * cword + (1 + cloze * cword | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 14 840.47 917.41 -406.24 812.47
## model 15 836.22 918.65 -403.11 806.22 6.2532 1 0.0124 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cloze:cword *** "
## [1] "remove . ~ . - cword "## refitting model(s) with ML (instead of REML)## Data: wordccompinput
## Models:
## model2: abssum ~ cloze + (1 + cloze * cword | sub)
## model: abssum ~ cloze + cword + (1 + cloze * cword | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 13 863.59 935.03 -418.79 837.59
## model 14 840.47 917.41 -406.24 812.47 25.113 1 5.407e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cword *** "
## [1] "remove . ~ . - cloze "## refitting model(s) with ML (instead of REML)## Data: wordccompinput
## Models:
## model2: abssum ~ (1 + cloze * cword | sub)
## model: abssum ~ cloze + (1 + cloze * cword | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 12 870.86 936.80 -423.43 846.86
## model 13 863.59 935.03 -418.79 837.59 9.2689 1 0.002331 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cloze *** "# create figure for paper
#mapping = aes_string(x="depth", y=dependMeasure, colour = "Condition",linetype="Condition")
drawERP(wordccompinput,dependMeasure,"Condition", span = 0.8, timeline.y=2.1)
if (!noslopes){ # when doing final version, change noslopes to FALSE
myggsave("img/cloze.png",width=figwidth,height=figheight) # control font size with height/width
}
n4P6cor = computeN4P6correlation(wordccompinput)
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.446316579469202"n4P6cor = computeN4P6correlation2(wordccompinput,"High Cloze","Low Cloze")
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.53014035585898"CUT AND PASTE INTO LATEX DOC click in box and use right arrow to see whole text
Cloze probabilities were manipulated in the language input on direct objects in active transitive items.
The same test word \textit{water} was used in each of the test sentences, but its cloze probability for each verb was manipulated. The same test word \textit{water} was used in each of the test sentences, but its cloze probability for each verb was manipulated. It occurred 60\% of the time with the verb \textit{drink} (High Cloze), 15\% with the verb \textit{taste} (Medium Cloze), and 4\% of the time with the verb \textit{take} (Low Cloze).
This is an example of the stimuli used in the model
\begin{tabular}{l|l}
\hline
Condition & Example\\
\hline
High Cloze & a teacher was drink -ing the water . .\\
\hline
Medium Cloze & a teacher was taste -ing the water . .\\
\hline
Low Cloze & a teacher was take -ing the water . .\\
\hline
\end{tabular}
The maximal model for the data had random slopes + (1 + cloze*cword |sub)
Mixed Model Results
a main effect of cloze ($\beta$=0.41, SE=0.058, $\chi^2$(1)=9.27, p=0.0023)
a main effect of cword ($\beta$=0.68, SE=0.059, $\chi^2$(1)=25.11, p<0.001)
@@an interaction of cloze:cword ($\beta$=0.19, SE=0.069, $\chi^2$(1)=6.25, p=0.012)
Federmeier et al. 2007
write("feder", stderr())
feder.df = subset(adultdata.df, str_detect(adultdata.df$Condition,"(EXPECT)") & wordcat %in% c('NOUNI'))
cond = str_split_fixed(feder.df$Condition,";",2)
feder.df$Constraint = cond[,1]
feder.df$Expectation = cond[,2]
feder.df$Constraint = factor(feder.df$Constraint,labels=c("Strong","Weak"))
feder.df$Expectation = factor(feder.df$Expectation,labels=c("Expected","Unexpected"))
nc = length(feder.df)
feder.df = feder.df[,c(nc-1,nc,1:(nc-2))] # move new columns to front
print(head(feder.df[feder.df$measure=="target",c("word","wordcat","Condition","Example","tick")])) # check## word wordcat Condition
## 114492 tea NOUNI STRONGCONS;EXPECTED
## 114596 wine NOUNI WEAKCONS;EXPECTED
## 114700 water NOUNI STRONGCONS;UNEXPECTED
## 114804 water NOUNI WEAKCONS;UNEXPECTED
## 115961 tea NOUNI STRONGCONS;EXPECTED
## 116091 wine NOUNI WEAKCONS;EXPECTED
## Example tick
## 114492 a wife will sip the tea . . 6
## 114596 a wife will sniff the wine . . 6
## 114700 a wife will sip the water . . 6
## 114804 a wife will sniff the water . . 6
## 115961 the brother -s were sip -ing the tea . . 8
## 116091 the brother -s were sniff -ing the wine . . 8showWordOutError("coffee,tea,wine,water,cake",c("Expectation","Constraint","measure"),feder.df)
computeClozeERPCorr(feder.df,dependMeasure,"word")## [1] "Correlation output prediction and abssum error of word is -0.999999999999786"corrHiddenOutputInputDeriv(feder.df,"Expectation",dim2="Constraint",layername="word")## Condition
## ACTCONT CATCONT CATVIOL
## 0 0 0
## CONG GARDAMB GARDSC
## 0 0 0
## HIGHCLOZE INCOH LOWCLOZE
## 0 0 0
## PASSCONT PLUR;CONT PLUR;VIOL
## 0 0 0
## SEMATTR SING;CONT SING;VIOL
## 0 0 0
## STRONG STRONGCONS;EXPECTED STRONGCONS;UNEXPECTED
## 0 8 8
## SUBCONT SUBVIOL TENSECONT
## 0 0 0
## TENSEVIOL TEST-PRED TEST-UNPRED
## 0 0 0
## WEAK WEAKCONS;EXPECTED WEAKCONS;UNEXPECTED
## 0 8 8
## ́ZEROCLOZE
## 0
# this is for creating ERP figure
nooutput = subset(feder.df, ! measure %in% layerNotERP )
end = which(names(nooutput)=="depth")-1
first = which(names(nooutput)=="none")
norawfeder.df = nooutput[,-(first:end)] # remove word specific activations/error
dvlen = length(names(norawfeder.df))- which(names(norawfeder.df)=="depth")
longnoraw= melt(norawfeder.df,id.var=1:(length(norawfeder.df)-dvlen)) # make long format
print(folder) # we use folder rename html file## [1] "sim80-s0-Sim116-alt10-50/deriv0.1vcoul1-100000"# add target layer for figure
#mapping = aes(x=depth, y=value, colour = Constraint,linetype=Expectation)
p1= drawERP(longnoraw,"value","Constraint","Expectation","measure","variable", timesize=-1)
p1 = p1 + facet_wrap( ~ measure+variable,scales = "free",nrow=1)
p1
# get example for table
tablesub = subset(norawfeder.df,layer == "word" & measure == mainMeasure)
exampleTable = tablesub[1:4,c("Constraint","Expectation","Example")]
tabletext = kable(exampleTable,format="latex",row.names=FALSE)
print(exampleTable)## Constraint Expectation Example
## 114494 Strong Expected a wife will sip the tea . .
## 114598 Weak Expected a wife will sniff the wine . .
## 114702 Strong Unexpected a wife will sip the water . .
## 114806 Weak Unexpected a wife will sniff the water . .interactionText=""
wordccompinput = subset(norawfeder.df,measure == mainMeasure & layer %in% c("word",p600layer) )
if (nlevels(adultdata.df$sub) > 4){ # do mixed with more than 9 subj
print("run models")
# our analysis uses this set
wordccompinput$cstrong = ifelse(wordccompinput$Constraint=="Strong",-0.5,0.5)
wordccompinput$cexpected = ifelse(wordccompinput$Expectation=="Expected",-0.5,0.5)
wordccompinput$cword = ifelse(wordccompinput$layer=="word",0.5,-0.5)
# - Constraint:Expectation:layer - Constraint:Expectation - Constraint:layer - Expectation:layer
# this is the main mixed model with centered variables
randomeff = " + (1 + Constraint*Expectation*layer - Constraint:Expectation:layer - Constraint:layer | sub)"
formu = as.formula(paste(dependMeasure,"~ cexpected + cstrong + cword + Constraint+Expectation+layer + sub" ))
if (noslopes){ randomeff = " + (1 | sub)"}
formu = as.formula(paste(dependMeasure,"~ cexpected*cstrong*cword",randomeff ))
omnimodel = lmer(formu, wordccompinput,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodel))
modellist = modelComparison(omnimodel) # create anova table by terms subtraction
randomeff = " + (1 + Constraint*Expectation | sub)"
wordccompinputword = subset(wordccompinput, layer == "word")
formu = as.formula(paste(dependMeasure,"~ cexpected*cstrong",randomeff ))
omnimodelword = lmer(formu, wordccompinputword,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodelword))
modellistword = modelComparison(omnimodelword) # create anova table by terms subtraction
wordccompinputhidden= subset(wordccompinput, layer == "hidden")
omnimodelhid = lmer(formu, wordccompinputhidden,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodelhid))
modellisthid = modelComparison(omnimodelhid) # create anova table by terms subtraction
# this is the equivalent model using factors that we need for posthocs
randomeff2 = " + (1 + Expectation*layer | sub)"
formu2 = as.formula(paste(dependMeasure,"~ Expectation*layer",randomeff2 ))
omnimodelFactor = lmer(formu2, wordccompinput)
# print(summary(omnimodelFactor))
model.lsmobj <- lsmeans(omnimodelFactor, ~ Expectation | layer)
posthocs = summary(as.glht(pairs(model.lsmobj)))
print(posthocs)
interactionText = printInteraction(posthocs)
randomeff2 = " + (1 + Constraint*Expectation*layer | sub)"
formu2 = as.formula(paste(dependMeasure,"~ Constraint*Expectation*layer",randomeff2 ))
omnimodelFactor = lmer(formu2, wordccompinput)
# print(summary(omnimodelFactor))
model.lsmobj <- lsmeans(omnimodelFactor, ~ Constraint*Expectation | layer)
posthocs = summary(as.glht(pairs(model.lsmobj)))
print(posthocs)
interactionText = printInteraction(posthocs)
}## [1] "run models"
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## abssum ~ cexpected * cstrong * cword + (1 + Constraint * Expectation *
## layer - Constraint:Expectation:layer - Constraint:layer | sub)
## Data: wordccompinput
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: 1024.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.2155 -0.4290 -0.0236 0.3064 6.7321
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.05946 0.2438
## ConstraintWeak 0.11114 0.3334 -0.93
## ExpectationUnexpected 0.12539 0.3541 -0.86
## layerhidden 0.03781 0.1944 -0.60
## ConstraintWeak:ExpectationUnexpected 0.34447 0.5869 0.69
## ExpectationUnexpected:layerhidden 0.27877 0.5280 0.13
## Residual 0.08333 0.2887
##
##
##
## 0.79
## 0.49 0.55
## -0.72 -0.93 -0.58
## -0.08 0.37 0.12 -0.60
##
## Number of obs: 2400, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.03553 0.04704 22.013
## cexpected 0.54516 0.08122 6.712
## cstrong 0.01401 0.07609 0.184
## cword 0.71520 0.11030 6.484
## cexpected:cstrong -0.75614 0.18709 -4.042
## cexpected:cword 0.34147 0.16862 2.025
## cstrong:cword 0.40188 0.02357 17.051
## cexpected:cstrong:cword 0.01286 0.04714 0.273
##
## Correlation of Fixed Effects:
## (Intr) cxpctd cstrng cword cxpctd:cs cxpctd:cw cstrn:
## cexpected 0.765
## cstrong -0.742 -0.662
## cword -0.935 -0.768 0.654
## cxpctd:cstr -0.699 -0.748 0.221 0.773
## cxpctd:cwrd -0.926 -0.841 0.838 0.818 0.592
## cstrng:cwrd 0.000 0.000 0.000 0.000 0.000 0.000
## cxpctd:cst: 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## [1] "remove . ~ . - cexpected:cstrong:cword "## refitting model(s) with ML (instead of REML)## Data: wordccompinput
## Models:
## model2: abssum ~ cexpected + cstrong + cword + (1 + Constraint * Expectation *
## model2: layer - Constraint:Expectation:layer - Constraint:layer |
## model2: sub) + cexpected:cstrong + cexpected:cword + cstrong:cword
## model: abssum ~ cexpected * cstrong * cword + (1 + Constraint * Expectation *
## model: layer - Constraint:Expectation:layer - Constraint:layer |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 29 1046.3 1214.0 -494.17 988.33
## model 30 1048.3 1221.8 -494.13 988.26 0.0745 1 0.7848
## [1] "########## Above comparison for cexpected:cstrong:cword "
## [1] "remove . ~ . - cstrong:cword "## refitting model(s) with ML (instead of REML)## Data: wordccompinput
## Models:
## model2: abssum ~ cexpected + cstrong + cword + (1 + Constraint * Expectation *
## model2: layer - Constraint:Expectation:layer - Constraint:layer |
## model2: sub) + cexpected:cstrong + cexpected:cword
## model: abssum ~ cexpected + cstrong + cword + (1 + Constraint * Expectation *
## model: layer - Constraint:Expectation:layer - Constraint:layer |
## model: sub) + cexpected:cstrong + cexpected:cword + cstrong:cword
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 28 1318.8 1480.7 -631.38 1262.75
## model 29 1046.3 1214.0 -494.17 988.33 274.42 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cstrong:cword *** "
## [1] "remove . ~ . - cexpected:cword "## refitting model(s) with ML (instead of REML)## Data: wordccompinput
## Models:
## model2: abssum ~ cexpected + cstrong + cword + (1 + Constraint * Expectation *
## model2: layer - Constraint:Expectation:layer - Constraint:layer |
## model2: sub) + cexpected:cstrong
## model: abssum ~ cexpected + cstrong + cword + (1 + Constraint * Expectation *
## model: layer - Constraint:Expectation:layer - Constraint:layer |
## model: sub) + cexpected:cstrong + cexpected:cword
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 27 1320.5 1476.7 -633.25 1266.5
## model 28 1318.8 1480.7 -631.38 1262.8 3.7491 1 0.05284 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cexpected:cword "
## [1] "remove . ~ . - cexpected:cstrong "## refitting model(s) with ML (instead of REML)## Data: wordccompinput
## Models:
## model2: abssum ~ cexpected + cstrong + cword + (1 + Constraint * Expectation *
## model2: layer - Constraint:Expectation:layer - Constraint:layer |
## model2: sub)
## model: abssum ~ cexpected + cstrong + cword + (1 + Constraint * Expectation *
## model: layer - Constraint:Expectation:layer - Constraint:layer |
## model: sub) + cexpected:cstrong
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 26 1333.0 1483.3 -640.48 1281.0
## model 27 1320.5 1476.7 -633.25 1266.5 14.46 1 0.0001432 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cexpected:cstrong *** "
## [1] "remove . ~ . - cword "## refitting model(s) with ML (instead of REML)## Data: wordccompinput
## Models:
## model2: abssum ~ cexpected + cstrong + (1 + Constraint * Expectation *
## model2: layer - Constraint:Expectation:layer - Constraint:layer |
## model2: sub)
## model: abssum ~ cexpected + cstrong + cword + (1 + Constraint * Expectation *
## model: layer - Constraint:Expectation:layer - Constraint:layer |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 25 1348.1 1492.7 -649.07 1298.1
## model 26 1333.0 1483.3 -640.48 1281.0 17.173 1 3.413e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cword *** "
## [1] "remove . ~ . - cstrong "## refitting model(s) with ML (instead of REML)## Data: wordccompinput
## Models:
## model2: abssum ~ cexpected + (1 + Constraint * Expectation * layer -
## model2: Constraint:Expectation:layer - Constraint:layer | sub)
## model: abssum ~ cexpected + cstrong + (1 + Constraint * Expectation *
## model: layer - Constraint:Expectation:layer - Constraint:layer |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 24 1353.6 1492.4 -652.80 1305.6
## model 25 1348.1 1492.7 -649.07 1298.1 7.4682 1 0.00628 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cstrong *** "
## [1] "remove . ~ . - cexpected "## refitting model(s) with ML (instead of REML)## Data: wordccompinput
## Models:
## model2: abssum ~ (1 + Constraint * Expectation * layer - Constraint:Expectation:layer -
## model2: Constraint:layer | sub)
## model: abssum ~ cexpected + (1 + Constraint * Expectation * layer -
## model: Constraint:Expectation:layer - Constraint:layer | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 23 1351.8 1484.8 -652.91 1305.8
## model 24 1353.6 1492.4 -652.80 1305.6 0.2126 1 0.6448
## [1] "########## Above comparison for cexpected "
## Linear mixed model fit by REML ['lmerMod']
## Formula: abssum ~ cexpected * cstrong + (1 + Constraint * Expectation |
## sub)
## Data: wordccompinputword
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: -3271.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1401 -0.3910 0.0104 0.3631 5.8163
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.052077 0.22820
## ConstraintWeak 0.088918 0.29819 -0.90
## ExpectationUnexpected 0.074695 0.27330 -0.98
## ConstraintWeak:ExpectationUnexpected 0.136510 0.36947 0.89
## Residual 0.003253 0.05704
##
##
##
## 0.88
## -0.99 -0.87
##
## Number of obs: 1200, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.39313 0.01974 70.576
## cexpected 0.71589 0.04619 15.499
## cstrong 0.21495 0.03761 5.715
## cexpected:cstrong -0.74970 0.11702 -6.406
##
## Correlation of Fixed Effects:
## (Intr) cxpctd cstrng
## cexpected -0.451
## cstrong -0.337 0.421
## cxpctd:cstr 0.482 -0.358 -0.926
## [1] "remove . ~ . - cexpected:cstrong "## refitting model(s) with ML (instead of REML)## Data: wordccompinputword
## Models:
## model2: abssum ~ cexpected + cstrong + (1 + Constraint * Expectation |
## model2: sub)
## model: abssum ~ cexpected * cstrong + (1 + Constraint * Expectation |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 14 -3246.4 -3175.1 1637.2 -3274.4
## model 15 -3261.5 -3185.2 1645.8 -3291.5 17.157 1 3.442e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cexpected:cstrong *** "
## [1] "remove . ~ . - cstrong "## refitting model(s) with ML (instead of REML)## Data: wordccompinputword
## Models:
## model2: abssum ~ cexpected + (1 + Constraint * Expectation | sub)
## model: abssum ~ cexpected + cstrong + (1 + Constraint * Expectation |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 13 -3248.3 -3182.1 1637.2 -3274.3
## model 14 -3246.4 -3175.1 1637.2 -3274.4 0.0656 1 0.7978
## [1] "########## Above comparison for cstrong "
## [1] "remove . ~ . - cexpected "## refitting model(s) with ML (instead of REML)## Data: wordccompinputword
## Models:
## model2: abssum ~ (1 + Constraint * Expectation | sub)
## model: abssum ~ cexpected + (1 + Constraint * Expectation | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 12 -3233.6 -3172.5 1628.8 -3257.6
## model 13 -3248.3 -3182.1 1637.2 -3274.3 16.738 1 4.29e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cexpected *** "
## Linear mixed model fit by REML ['lmerMod']
## Formula: abssum ~ cexpected * cstrong + (1 + Constraint * Expectation |
## sub)
## Data: wordccompinputhidden
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: 1037.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.4447 -0.4338 -0.1283 0.3220 4.3005
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.05572 0.2361
## ConstraintWeak 0.17206 0.4148 -0.77
## ExpectationUnexpected 0.77526 0.8805 0.14
## ConstraintWeak:ExpectationUnexpected 0.73853 0.8594 0.07
## Residual 0.12530 0.3540
##
##
##
## 0.07
## -0.38 -0.94
##
## Number of obs: 1200, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.6779 0.1005 6.744
## cexpected 0.3744 0.1588 2.358
## cstrong -0.1869 0.1497 -1.249
## cexpected:cstrong -0.7626 0.2748 -2.775
##
## Correlation of Fixed Effects:
## (Intr) cxpctd cstrng
## cexpected 0.867
## cstrong -0.818 -0.900
## cxpctd:cstr -0.816 -0.785 0.565
## [1] "remove . ~ . - cexpected:cstrong "## refitting model(s) with ML (instead of REML)## Data: wordccompinputhidden
## Models:
## model2: abssum ~ cexpected + cstrong + (1 + Constraint * Expectation |
## model2: sub)
## model: abssum ~ cexpected * cstrong + (1 + Constraint * Expectation |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 14 1058.9 1130.2 -515.45 1030.9
## model 15 1054.7 1131.1 -512.36 1024.7 6.1817 1 0.01291 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cexpected:cstrong *** "
## [1] "remove . ~ . - cstrong "## refitting model(s) with ML (instead of REML)## Data: wordccompinputhidden
## Models:
## model2: abssum ~ cexpected + (1 + Constraint * Expectation | sub)
## model: abssum ~ cexpected + cstrong + (1 + Constraint * Expectation |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 13 1057.0 1123.2 -515.50 1031.0
## model 14 1058.9 1130.2 -515.45 1030.9 0.0893 1 0.7651
## [1] "########## Above comparison for cstrong "
## [1] "remove . ~ . - cexpected "## refitting model(s) with ML (instead of REML)## Data: wordccompinputhidden
## Models:
## model2: abssum ~ (1 + Constraint * Expectation | sub)
## model: abssum ~ cexpected + (1 + Constraint * Expectation | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 12 1056.1 1117.2 -516.06 1032.1
## model 13 1057.0 1123.2 -515.50 1031.0 1.1169 1 0.2906
## [1] "########## Above comparison for cexpected "
## $`layer = word`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Expected - Unexpected == 0 -0.71589 0.04479 -15.98 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
##
## $`layer = hidden`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Expected - Unexpected == 0 -0.3744 0.1588 -2.358 0.0184 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
##
## $`layer = word`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value
## Strong,Expected - Weak,Expected == 0 -0.58981 0.09455 -6.238
## Strong,Expected - Strong,Unexpected == 0 -1.09075 0.08658 -12.598
## Strong,Expected - Weak,Unexpected == 0 -0.93085 0.07160 -13.002
## Weak,Expected - Strong,Unexpected == 0 -0.50094 0.04707 -10.643
## Weak,Expected - Weak,Unexpected == 0 -0.34104 0.06058 -5.630
## Strong,Unexpected - Weak,Unexpected == 0 0.15990 0.03107 5.146
## Pr(>|z|)
## Strong,Expected - Weak,Expected == 0 < 1e-06 ***
## Strong,Expected - Strong,Unexpected == 0 < 1e-06 ***
## Strong,Expected - Weak,Unexpected == 0 < 1e-06 ***
## Weak,Expected - Strong,Unexpected == 0 < 1e-06 ***
## Weak,Expected - Weak,Unexpected == 0 < 1e-06 ***
## Strong,Unexpected - Weak,Unexpected == 0 1.07e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
##
## $`layer = hidden`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value
## Strong,Expected - Weak,Expected == 0 -0.194357 0.134313 -1.447
## Strong,Expected - Strong,Unexpected == 0 -0.755710 0.278231 -2.716
## Strong,Expected - Weak,Unexpected == 0 -0.187500 0.069325 -2.705
## Weak,Expected - Strong,Unexpected == 0 -0.561353 0.298719 -1.879
## Weak,Expected - Weak,Unexpected == 0 0.006857 0.098873 0.069
## Strong,Unexpected - Weak,Unexpected == 0 0.568210 0.253095 2.245
## Pr(>|z|)
## Strong,Expected - Weak,Expected == 0 0.3964
## Strong,Expected - Strong,Unexpected == 0 0.0249 *
## Strong,Expected - Weak,Unexpected == 0 0.0256 *
## Weak,Expected - Strong,Unexpected == 0 0.1878
## Weak,Expected - Weak,Unexpected == 0 0.9998
## Strong,Unexpected - Weak,Unexpected == 0 0.0853 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)# create figure for paper
wordccompinput$layer= factor(wordccompinput$layer)
#mapping = aes_string(x="depth", y=dependMeasure, colour = "Constraint",linetype="Expectation")
drawERP(wordccompinput,dependMeasure,"Expectation","Constraint" , timeline.y=2)
if (!noslopes){ # when doing final version, change noslopes to FALSE
myggsave("img/fed2007.png",width=figwidth,height=figheight) # control font size with height/width
}
n4P6cor = computeN4P6correlation(wordccompinput)
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.350002016138098"CUT AND PASTE INTO LATEX DOC click in box and use right arrow to see whole text
Sentential constraint was varied by manipulating the verb. The strong sentential constraint verb \textit{sip} occurred with \textit{tea} 61\% of the time and \textit{water} 10\% of the time. The weak sentential constraint verb \textit{sniff} occurred with \textit{wine} and with \textit{water} 15\% of the time. Word expectation can be examined by comparing the unexpected word \textit{water} and the expected words \textit{tea} and \textit{wine}.
This is an example of the stimuli used in the model
\begin{tabular}{l|l|l}
\hline
Constraint & Expectation & Example\\
\hline
Strong & Expected & a wife will sip the tea . .\\
\hline
Weak & Expected & a wife will sniff the wine . .\\
\hline
Strong & Unexpected & a wife will sip the water . .\\
\hline
Weak & Unexpected & a wife will sniff the water . .\\
\hline
\end{tabular}
The maximal model for the data had random slopes + (1 + Constraint*Expectation | sub)
Mixed Model Results
no main effect of cexpected ($\beta$=0.55, SE=0.081, $\chi^2$(1)=0.21, p=0.645)
a main effect of cstrong ($\beta$=0.014, SE=0.076, $\chi^2$(1)=7.47, p=0.0063)
a main effect of cword ($\beta$=0.72, SE=0.11, $\chi^2$(1)=17.17, p<0.001)
an interaction of cexpected:cstrong ($\beta$=-0.76, SE=0.19, $\chi^2$(1)=14.46, p<0.001)
no interaction of cexpected:cword ($\beta$=0.34, SE=0.17, $\chi^2$(1)=3.75, p=0.053)
an interaction of cstrong:cword ($\beta$=0.4, SE=0.024, $\chi^2$(1)=274.42, p<0.001)
@@no interaction of cexpected:cstrong:cword ($\beta$=0.013, SE=0.047, $\chi^2$(1)=0.07, p=0.785)
Mixed Model Results
a main effect of cexpected ($\beta$=0.72, SE=0.046, $\chi^2$(1)=16.74, p<0.001)
no main effect of cstrong ($\beta$=0.21, SE=0.038, $\chi^2$(1)=0.07, p=0.798)
@@an interaction of cexpected:cstrong ($\beta$=-0.75, SE=0.12, $\chi^2$(1)=17.16, p<0.001)
Mixed Model Results
no main effect of cexpected ($\beta$=0.37, SE=0.16, $\chi^2$(1)=1.12, p=0.291)
no main effect of cstrong ($\beta$=-0.19, SE=0.15, $\chi^2$(1)=0.09, p=0.765)
@@an interaction of cexpected:cstrong ($\beta$=-0.76, SE=0.27, $\chi^2$(1)=6.18, p=0.013)
Posthoc tests for each layer
@@Posthoc: There was a significant difference for layer = word, diff = 0.5898, t(0)=6.24, p<0.001
@@Posthoc: There was no difference for layer = hidden, p=0.3964
ratio layer = hidden 0.329603255340793 rev 3.03395061728395
Posthoc tests for each layer
@@Posthoc: There was a significant difference for layer = word, diff = 0.5898, t(0)=6.24, p<0.001
@@Posthoc: There was no difference for layer = hidden, p=0.0249
ratio layer = hidden 0.329603255340793 rev 3.03395061728395
Posthoc tests for each layer
@@Posthoc: There was a significant difference for layer = word, diff = 0.5898, t(0)=6.24, p<0.001
@@Posthoc: There was no difference for layer = hidden, p=0.0256
ratio layer = hidden 0.329603255340793 rev 3.03395061728395
Posthoc tests for each layer
@@Posthoc: There was a significant difference for layer = word, diff = 0.5898, t(0)=6.24, p<0.001
@@Posthoc: There was no difference for layer = hidden, p=0.1878
ratio layer = hidden 0.329603255340793 rev 3.03395061728395
Posthoc tests for each layer
@@Posthoc: There was a significant difference for layer = word, diff = 0.5898, t(0)=6.24, p<0.001
@@Posthoc: There was no difference for layer = hidden, p=0.9998
ratio layer = hidden 0.329603255340793 rev 3.03395061728395
Posthoc tests for each layer
@@Posthoc: There was a significant difference for layer = word, diff = 0.5898, t(0)=6.24, p<0.001
@@Posthoc: There was no difference for layer = hidden, p=0.0853
ratio layer = hidden 0.329603255340793 rev 3.03395061728395
Position Van Petten and Kutas
write("vanpetten", stderr())
pos.df = subset(adultdata.df, str_detect(adultdata.df$Condition,"(CONG|INCOH)")) # & tick < 10)
print(head(pos.df[pos.df$measure=="target",c("word","wordcat","Condition","Example","tick")])) # check## word wordcat Condition Example
## 106913 man NOUNA CONG man -s throw -ed the toy to the driver . .
## 106926 -s PL CONG man -s throw -ed the toy to the driver . .
## 106939 throw VERBD CONG man -s throw -ed the toy to the driver . .
## 106952 -ed ED CONG man -s throw -ed the toy to the driver . .
## 106965 the DET CONG man -s throw -ed the toy to the driver . .
## 106978 toy NOUNI CONG man -s throw -ed the toy to the driver . .
## tick
## 106913 1
## 106926 2
## 106939 3
## 106952 4
## 106965 5
## 106978 6postargword.df = subset(pos.df, layer %in% c("word") & measure %in% c("output",mainMeasure) & contfunc=="C" & wordcat != "VERBD" & Condition == "CONG")
postargword.df$verb = str_match(postargword.df$Example,"(give|throw|send|lend)")[,2]
meanbyword = aggregate(cbind(man,kite,cake) ~ tick + measure + verb + Condition,postargword.df, mean)
meanbywordlong= melt(meanbyword,id.var=1:4)
p1 = ggplot(meanbywordlong, aes(x=verb, y=value, fill = variable))
p1 = p1 + geom_bar(stat="identity",position="dodge")
#p1 = p1 + facet_wrap(~ measure + tick,scales = "free",ncol=4)
p1 = p1 + facet_grid( tick ~ measure + Condition,scales = "free_y")
p1 = p1 + theme(legend.position="bottom")
p1
# this is for creating ERP figure
nooutput = subset(pos.df, ! measure %in% layerNotERP)
end = which(names(nooutput)=="depth")-1
first = which(names(nooutput)=="none")
norawpos.df = nooutput[,-(first:end)] # remove word specific activations/error
norawpos.df$Condition = factor(norawpos.df$Condition,labels=c("Congruent","Syntactic"))
dvlen = length(names(norawpos.df))- which(names(norawpos.df)=="depth")
longnoraw= melt(norawpos.df,id.var=1:(length(norawpos.df)-dvlen)) # make long format
# add target layer for figure
longnoraw2 = subset(longnoraw, measure == mainMeasure & contfunc == "C") #& wordcat != "VERBD"
#longnoraw2$contfunc = paste(longnoraw2$tick,longnoraw2$contfunc)
#mapping = aes(x=depth, y=value, colour = Condition,linetype=Condition)
p1= drawERP(longnoraw2,"value","Condition","variable","tick", timesize=-1)
p1 = p1 + facet_grid(tick ~ variable ,scales = "free")
#p1 = p1 + facet_wrap( ~ tick+variable)
p1
interactionText=""
# get example for table
tablesub = subset(norawpos.df,layer == "word" & measure == mainMeasure & tick == 2)
exampleTable = tablesub[1:4,c("Condition","Example")]
tabletext = kable(exampleTable,format="latex",row.names=FALSE)
print(exampleTable)## Condition Example
## 106928 Congruent man -s throw -ed the toy to the driver . .
## 107071 Syntactic steak -s throw -ed the man to the kite . .
## 107630 Congruent the father was lend -ing a sister the kite -s . .
## 107786 Syntactic the ball was lend -ing a coffee the woman -s . .posinput = subset(norawpos.df,measure == mainMeasure & layer %in% c("word",p600layer) & word != "." & contfunc=="C" & wordcat %in% c("NOUNA","NOUNI"))
#posinput = subset(posinput, !str_detect(Example,"-par") )
if (nlevels(adultdata.df$sub) > 4){ # do mixed with more than 9 subj
print("run models")
# our analysis uses this set
posinput$dcoher = ifelse(posinput$Condition=="Congruent",1,0)
posinput$ccoher = ifelse(posinput$Condition=="Congruent",0.5,-0.5)
posinput$cword = ifelse(posinput$layer=="word",0.5,-0.5)
posinput$ctick = posinput$tick - mean(posinput$tick)
posinput$ccontent = ifelse(posinput$contfunc=="C",0.5,-0.5)
#& wordcat != "VERBD"
posinputnoun = subset(posinput, layer == "word")
randomeff = paste(" + (1 + dcoher*ctick |sub)")
if (noslopes){ randomeff = " + (1 | sub)"}
formu = as.formula(paste(dependMeasure,"~ dcoher*ctick",randomeff ))
omnimodel = lmer(formu, posinputnoun,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodel))
modellist = modelComparison(omnimodel) # create anova table by terms subtraction
# create figure for paper
}## [1] "run models"
## Linear mixed model fit by REML ['lmerMod']
## Formula: abssum ~ dcoher * ctick + (1 + dcoher * ctick | sub)
## Data: posinputnoun
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: -7090.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.5023 -0.4530 0.0886 0.4666 3.6429
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 3.644e-06 0.001909
## dcoher 2.932e-05 0.005414 -1.00
## ctick 7.671e-08 0.000277 1.00 -1.00
## dcoher:ctick 1.797e-06 0.001341 -1.00 1.00 -1.00
## Residual 1.103e-03 0.033215
## Number of obs: 1800, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.9710800 0.0012610 1563.05
## dcoher -0.0583900 0.0023202 -25.17
## ctick 0.0040394 0.0003498 11.55
## dcoher:ctick -0.0109311 0.0006396 -17.09
##
## Correlation of Fixed Effects:
## (Intr) dcoher ctick
## dcoher -0.772
## ctick 0.120 -0.185
## dcoher:ctck -0.317 0.489 -0.679
## [1] "remove . ~ . - dcoher:ctick "## refitting model(s) with ML (instead of REML)## Data: posinputnoun
## Models:
## model2: abssum ~ dcoher + ctick + (1 + dcoher * ctick | sub)
## model: abssum ~ dcoher * ctick + (1 + dcoher * ctick | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 14 -7077.1 -7000.1 3552.5 -7105.1
## model 15 -7110.9 -7028.5 3570.4 -7140.9 35.824 1 2.16e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for dcoher:ctick *** "
## [1] "remove . ~ . - ctick "## refitting model(s) with ML (instead of REML)## Data: posinputnoun
## Models:
## model2: abssum ~ dcoher + (1 + dcoher * ctick | sub)
## model: abssum ~ dcoher + ctick + (1 + dcoher * ctick | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 13 -7079.0 -7007.5 3552.5 -7105.0
## model 14 -7077.1 -7000.1 3552.5 -7105.1 0.0863 1 0.769
## [1] "########## Above comparison for ctick "
## [1] "remove . ~ . - dcoher "## refitting model(s) with ML (instead of REML)## Data: posinputnoun
## Models:
## model2: abssum ~ (1 + dcoher * ctick | sub)
## model: abssum ~ dcoher + (1 + dcoher * ctick | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 12 -7062.3 -6996.4 3543.2 -7086.3
## model 13 -7079.0 -7007.5 3552.5 -7105.0 18.673 1 1.552e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for dcoher *** "# & wordcat != "VERBD"
longnoraw2 = subset(posinput, layer == "word" )
form = as.formula(paste(dependMeasure," ~ tick + Condition",sep=""))
mdf = aggregate(form, longnoraw2, mean)
mapping = aes_string(x="tick", y=dependMeasure, colour = "Condition",linetype="Condition")
p1 = ggplot(mdf,mapping=mapping)
#p1 = p1 + geom_line() + geom_point()
p1 = p1 + geom_point()
p1 = p1 + stat_smooth(method="lm",se=F)
p1 = p1 + scale_colour_grey()
# p1 = p1 + scale_colour_brewer( palette="Set1")
p1 = p1 + scale_x_continuous(breaks=1:12)
p1 = p1 + theme_bw()
p1 = p1 +theme(legend.background = element_rect(size=.4, color="grey80"))
p1=p1 + theme(legend.position="bottom")+xlab("Sentence Position")+ylab("Sum Abs. Error")
p1 = p1 + theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"))
cong = str_split(exampleTable$Example[1]," ")[[1]]
inco = str_split(exampleTable$Example[2]," ")[[1]]
coindf = data.frame(tick=rep(1:9,2),word=c(cong[2:10],inco[2:10]),Condition=rep(c("Congruent","Incoherent"),each=9),dependMeasure=rep(c(1.9,1.97),each=9))
# p1=p1+geom_text(aes(x=tick, y= dependMeasure,label=word,colour = Condition),coindf,size=3)
p1
if (!noslopes){ # when doing final version, change noslopes to FALSE
myggsave("img/vanpetten.png",width=figwidth,height=figheight) # control font size with height/width
}
n4P6cor = computeN4P6correlation(posinput)
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.434689569660451"CUT AND PASTE INTO LATEX DOC click in box and use right arrow to see whole text
Need stats to explain this
This is an example of the stimuli used in the model
\begin{tabular}{l|l}
\hline
Condition & Example\\
\hline
Congruent & man -s throw -ed the toy to the driver . .\\
\hline
Syntactic & steak -s throw -ed the man to the kite . .\\
\hline
Congruent & the father was lend -ing a sister the kite -s . .\\
\hline
Syntactic & the ball was lend -ing a coffee the woman -s . .\\
\hline
\end{tabular}
The maximal model for the data had random slopes + (1 + dcoher*ctick |sub)
Mixed Model Results
a main effect of dcoher ($\beta$=-0.058, SE=0.0023, $\chi^2$(1)=18.67, p<0.001)
no main effect of ctick ($\beta$=0.004, SE=0.00035, $\chi^2$(1)=0.09, p=0.769)
@@ @@an interaction of dcoher:ctick ($\beta$=-0.011, SE=0.00064, $\chi^2$(1)=35.82, p<0.001)
AGREEMENT Hagoort et al 1993
# This is Hagoort et al 1993
write("agree", stderr())
agree.df = subset(adultdata.df, str_detect(Condition,"(PLUR|SING)") )
verbpos = which(agree.df$wordcat == "VERBT")
#verb.df = agree.df[verbpos,]
wordcatsamelen = which(agree.df$wordcat!=agree.df$wordcat[1])[1]-1 # same number of rows for each category
postverb.df = agree.df[verbpos + wordcatsamelen,] # position after verb
postverb.df=postverb.df[!is.na(postverb.df$Condition),]
cond = str_split_fixed(postverb.df$Condition,";",2)
postverb.df$Number = cond[,1]
postverb.df$Agreement = cond[,2]
postverb.df$Number = factor(postverb.df$Number,labels=c("Plural","Singular"))
postverb.df$Agreement = factor(postverb.df$Agreement,labels=c("Control","Violation"))
nc = length(postverb.df)
postverb.df = postverb.df[,c(nc-1,nc,1:(nc-2))]
xtabs(~ Agreement + Number ,postverb.df)## Number
## Agreement Plural Singular
## Control 3000 3000
## Violation 3000 3000print(head(postverb.df[postverb.df$measure=="target",c("word","wordcat","Condition","Example","tick")])) # check## word wordcat Condition Example tick
## 107251 a DET PLUR;CONT the boy -s take a stick . . 5
## 107355 -ss SS PLUR;VIOL the boy -s take -ss a stick . . 5
## 107459 -ss SS SING;CONT the boy take -ss a stick . . 4
## 107563 a DET SING;VIOL the boy take a stick . . 4
## 110319 wine NOUNI PLUR;CONT the father -s taste wine -s . . 5
## 110423 -ss SS PLUR;VIOL the father -s taste -ss wine -s . . 5showWordOutError("the,X.ss,X.ing,X.ed,it",c("Agreement","measure","Condition"),postverb.df)
computeClozeERPCorr(postverb.df,dependMeasure,"hidden")## [1] "Correlation output prediction and abssum error of hidden is -0.706818447847079"corrHiddenOutputInputDeriv(postverb.df,"Agreement","Number")## Condition
## ACTCONT CATCONT CATVIOL
## 0 0 0
## CONG GARDAMB GARDSC
## 0 0 0
## HIGHCLOZE INCOH LOWCLOZE
## 0 0 0
## PASSCONT PLUR;CONT PLUR;VIOL
## 0 8 8
## SEMATTR SING;CONT SING;VIOL
## 0 8 8
## STRONG STRONGCONS;EXPECTED STRONGCONS;UNEXPECTED
## 0 0 0
## SUBCONT SUBVIOL TENSECONT
## 0 0 0
## TENSEVIOL TEST-PRED TEST-UNPRED
## 0 0 0
## WEAK WEAKCONS;EXPECTED WEAKCONS;UNEXPECTED
## 0 0 0
## ́ZEROCLOZE
## 0
# make erp figure with different measures
#mapping = aes(x=depth, y=value, colour = Number,linetype=Agreement)
#norawpostverbnotargout.df=removeWordsDrawERP(postverb.df,"value","Number","Agreement","variable",showFig=T)
norawpostverbnotargout.df=removeWordsDrawERP(postverb.df,"value","Agreement","Number","measure",iv5="variable",showFig=T)## Number Agreement sub epoch word wordcat contfunc
## 1 Plural Control s0 100000 a DET F
## 2 Plural Control s0 100000 a DET F
## 3 Plural Control s0 100000 a DET F
## 4 Plural Control s0 100000 a DET F
## 5 Plural Control s0 100000 a DET F
## 6 Plural Control s0 100000 a DET F
## Example Condition tick layer measure depth
## 1 the boy -s take a stick . . PLUR;CONT 5 word outderiv 1
## 2 the boy -s take a stick . . PLUR;CONT 5 word inputderiv 1
## 3 the boy -s take a stick . . PLUR;CONT 5 compress outderiv 2
## 4 the boy -s take a stick . . PLUR;CONT 5 compress inputderiv 2
## 5 the boy -s take a stick . . PLUR;CONT 5 hidden outderiv 3
## 6 the boy -s take a stick . . PLUR;CONT 5 hidden inputderiv 3
## variable value
## 1 abssum 6.151653
## 2 abssum 1.674885
## 3 abssum 49.612011
## 4 abssum 1.590005
## 5 abssum 38.645623
## 6 abssum 2.294385
interactionText="NO TEXT"
modellist2 = list()
difflex = 0
# get example for table
agreeinputword = subset(norawpostverbnotargout.df,measure == mainMeasure & layer == "word")
exampleTable = agreeinputword[1:4,c("Number","Agreement","Example")]
tabletext = kable(exampleTable,format="latex",row.names=FALSE)
print(exampleTable)## Number Agreement Example
## 107253 Plural Control the boy -s take a stick . .
## 107357 Plural Violation the boy -s take -ss a stick . .
## 107461 Singular Control the boy take -ss a stick . .
## 107565 Singular Violation the boy take a stick . .agreeinput = subset(norawpostverbnotargout.df,measure == mainMeasure & layer %in% c("word",p600layer) )
if (nlevels(adultdata.df$sub) > 4){ # do mixed with more than 9 subj
agreeinput$cplural = ifelse(agreeinput$Number=="Plural",0.5,-0.5)
agreeinput$cviolation = ifelse(agreeinput$Agreement=="Violation",0.5,-0.5)
agreeinput$cp6layer = ifelse(agreeinput$layer==p600layer,0.5,-0.5)
# this is the main mixed model with centered variables
randomeff = paste(" + (1 + cviolation*cp6layer |sub)") #Agreement * layer
if (noslopes){ randomeff = " + (1 | sub)"}
formu = as.formula(paste(dependMeasure,"~ cviolation*cp6layer",randomeff))
omnimodel = lmer(formu, agreeinput,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodel))
modellist2 = modelComparison(omnimodel) # create anova table
# this model is for posthocs
formu2 = as.formula(paste(dependMeasure,"~ Agreement*layer",randomeff))
omnimodelFactor = lmer(formu2, agreeinput)
#print(summary(omnimodelFactor))
model.lsmobj <- lsmeans(omnimodelFactor, ~ Agreement | layer)
posthocs = summary(as.glht(pairs(model.lsmobj)))
print(posthocs)
interactionText = printInteraction(posthocs)
difflex = abs(posthocs$`layer = word`$test$coefficients[[1]])
}## Linear mixed model fit by REML ['lmerMod']
## Formula: abssum ~ cviolation * cp6layer + (1 + cviolation * cp6layer |
## sub)
## Data: agreeinput
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: 7649.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4173 -0.3718 0.0001 0.2489 5.2085
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.1700 0.4123
## cviolation 0.6324 0.7953 0.97
## cp6layer 0.6334 0.7959 1.00 0.98
## cviolation:cp6layer 2.6656 1.6327 0.98 1.00 0.99
## Residual 1.3812 1.1752
## Number of obs: 2400, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.9028 0.1326 21.896
## cviolation 3.4526 0.2560 13.486
## cp6layer 2.3568 0.2562 9.198
## cviolation:cp6layer 5.8037 0.5251 11.052
##
## Correlation of Fixed Effects:
## (Intr) cviltn cp6lyr
## cviolation 0.934
## cp6layer 0.965 0.942
## cvltn:cp6ly 0.946 0.965 0.952
## [1] "remove . ~ . - cviolation:cp6layer "## refitting model(s) with ML (instead of REML)## Data: agreeinput
## Models:
## model2: abssum ~ cviolation + cp6layer + (1 + cviolation * cp6layer |
## model2: sub)
## model: abssum ~ cviolation * cp6layer + (1 + cviolation * cp6layer |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 14 7692.5 7773.4 -3832.2 7664.5
## model 15 7667.7 7754.4 -3818.8 7637.7 26.775 1 2.286e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation:cp6layer *** "
## [1] "remove . ~ . - cp6layer "## refitting model(s) with ML (instead of REML)## Data: agreeinput
## Models:
## model2: abssum ~ cviolation + (1 + cviolation * cp6layer | sub)
## model: abssum ~ cviolation + cp6layer + (1 + cviolation * cp6layer |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 13 7691.8 7767.0 -3832.9 7665.8
## model 14 7692.5 7773.4 -3832.2 7664.5 1.3736 1 0.2412
## [1] "########## Above comparison for cp6layer "
## [1] "remove . ~ . - cviolation "## refitting model(s) with ML (instead of REML)## Data: agreeinput
## Models:
## model2: abssum ~ (1 + cviolation * cp6layer | sub)
## model: abssum ~ cviolation + (1 + cviolation * cp6layer | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 12 7697.5 7766.9 -3836.7 7673.5
## model 13 7691.8 7767.0 -3832.9 7665.8 7.6503 1 0.005676 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation *** "
## $`layer = word`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Control - Violation == 0 -0.55076 0.06926 -7.953 1.78e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
##
## $`layer = hidden`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Control - Violation == 0 -6.354 0.514 -12.36 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method) # figure for paper
# mapping = aes_string(x="depth", y=dependMeasure, colour = "Agreement",linetype="Agreement")
drawERP(agreeinput,dependMeasure,"Agreement", span=1, timeline.y=5.2)
if (!noslopes){ # when doing final version, change noslopes to FALSE
myggsave("img/agree.png",width=figwidth,height=figheight)
}
n4P6cor = computeN4P6correlation(agreeinput)
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.706818447469042"n4P6cor = computeN4P6correlation2(agreeinput,"Control","Violation",cond="Agreement")
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.61678356390425"CUT AND PASTE INTO LATEX DOC click in box and use right arrow to see whole text
The model received active transitive input with singular subjects 75\% of the time and plural subjects 25\% of the time.
This is an example of the stimuli used in the model.
\begin{tabular}{l|l|l}
\hline
Number & Agreement & Example\\
\hline
Plural & Control & the boy -s take a stick . .\\
\hline
Plural & Violation & the boy -s take -ss a stick . .\\
\hline
Singular & Control & the boy take -ss a stick . .\\
\hline
Singular & Violation & the boy take a stick . .\\
\hline
\end{tabular}
The maximal model for the data had random slopes + (1 + cviolation*cp6layer |sub)
Mixed Model Results
a main effect of cviolation ($\beta$=3.5, SE=0.26, $\chi^2$(1)=7.65, p=0.0057)
no main effect of cp6layer ($\beta$=2.4, SE=0.26, $\chi^2$(1)=1.37, p=0.241)
@@an interaction of cviolation:cp6layer ($\beta$=5.8, SE=0.53, $\chi^2$(1)=26.78, p<0.001)
Posthoc tests for each layer
@@Posthoc: There was a significant difference for layer = word, diff = 0.5508, t(0)=7.95, p<0.001
@@Posthoc: There was a significant difference for layer = hidden, diff = 6.3545, t(0)=12.36, p<0.001
ratio layer = hidden 11.5368554829339 rev 0.0866787316075222Tense Allen et al. (2003)
write("tense", stderr())
tense.df = subset(adultdata.df, str_detect(Condition,"(TENSECONT|TENSEVIOL)") )
tense.df$Condition = factor(tense.df$Condition,labels=c("Control","Violation"))
verbpos = which(tense.df$wordcat == "VERBT")
wordcatsamelen = which(tense.df$wordcat!=tense.df$wordcat[1])[1]-1
verb.df = tense.df[verbpos,] # verb
postverb.df = tense.df[verbpos + wordcatsamelen,] # position after verb
postverb.df=postverb.df[!is.na(postverb.df$Condition),]
print(head(postverb.df[postverb.df$measure=="target",c("word","wordcat","Condition","Example","tick")]))## word wordcat Condition
## 134044 the DET Control
## 134161 -ed ED Violation
## 142130 the DET Control
## 142247 -ed ED Violation
## 150424 the DET Control
## 150554 -ed ED Violation
## Example tick
## 134044 a father will drink the beer -s . . 5
## 134161 a father will drink -ed the beer -s . . 5
## 142130 the driver will taste the wine -s . . 5
## 142247 the driver will taste -ed the wine -s . . 5
## 150424 the driver -s will taste the beer -s . . 6
## 150554 the driver -s will taste -ed the beer -s . . 6showWordOutError("the,a,X.ss,X.ing,X.ed,X.par",c("Condition","measure"),postverb.df)
corrHiddenOutputInputDeriv(postverb.df,"Condition")## Condition
## Control Violation
## 8 8
myggsave("img/actderivtense.png",width=6,height=6)# make erp figure with different measures
norawpostverbnotargout.df=removeWordsDrawERP(postverb.df,"value","Condition")## sub epoch word wordcat contfunc Example
## 1 s0 100000 the DET F a father will drink the beer -s . .
## 2 s0 100000 the DET F a father will drink the beer -s . .
## 3 s0 100000 the DET F a father will drink the beer -s . .
## 4 s0 100000 the DET F a father will drink the beer -s . .
## 5 s0 100000 the DET F a father will drink the beer -s . .
## 6 s0 100000 the DET F a father will drink the beer -s . .
## Condition tick layer measure depth variable value
## 1 Control 5 word outderiv 1 abssum 1.1872990
## 2 Control 5 word inputderiv 1 abssum 0.3155038
## 3 Control 5 compress outderiv 2 abssum 3.9494715
## 4 Control 5 compress inputderiv 2 abssum 0.1043950
## 5 Control 5 hidden outderiv 3 abssum 2.1760727
## 6 Control 5 hidden inputderiv 3 abssum 0.2120343interactionText="NO TEXT"
modellist2 = list()
# get example for table
tenseexampleword = subset(norawpostverbnotargout.df,measure == mainMeasure & layer == "word")
exampleTable = tenseexampleword[1:2,c("Condition","Example")]
tabletext = kable(exampleTable,format="latex",row.names=FALSE)
print(exampleTable)## Condition Example
## 134046 Control a father will drink the beer -s . .
## 134163 Violation a father will drink -ed the beer -s . . tenseinput = subset(norawpostverbnotargout.df,measure == mainMeasure & layer %in% c("word",p600layer) )
if (nlevels(adultdata.df$sub) > 4){ # do mixed with more than 9 subj
tenseinput$cviolation = ifelse(tenseinput$Condition=="Violation",0.5,-0.5) # one against two
tenseinput$cp600layer = ifelse(tenseinput$layer==p600layer,0.5,-0.5)
# this is the main mixed model with centered variables
randomeff = paste(" + (1 + cviolation*cp600layer |sub)") # Condition*layer
if (noslopes){ randomeff = " + (1 | sub)"}
formu = as.formula(paste(dependMeasure,"~ cviolation*cp600layer",randomeff))
omnimodel = lmer(formu, tenseinput,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodel))
modellist2 = modelComparison(omnimodel) # create anova table
# this model is for posthocs
formu2 = as.formula(paste(dependMeasure,"~ Condition*layer",randomeff))
omnimodelFactor = lmer(formu2, tenseinput)
#print(summary(omnimodelFactor))
model.lsmobj <- lsmeans(omnimodelFactor, ~ Condition | layer)
posthocs = summary(as.glht(pairs(model.lsmobj)))
print(posthocs)
interactionText = printInteraction(posthocs)
difflex = difflex + abs(posthocs$`layer = word`$test$coefficients[[1]])
}## Linear mixed model fit by REML ['lmerMod']
## Formula:
## abssum ~ cviolation * cp600layer + (1 + cviolation * cp600layer |
## sub)
## Data: tenseinput
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: 3647.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0214 -0.6388 0.0003 0.5295 7.0231
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.1447 0.3804
## cviolation 0.6981 0.8356 0.87
## cp600layer 0.5832 0.7637 0.98 0.95
## cviolation:cp600layer 2.4595 1.5683 0.94 0.99 0.99
## Residual 1.1648 1.0793
## Number of obs: 1200, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.7546 0.1243 22.170
## cviolation 3.1310 0.2715 11.534
## cp600layer 2.3415 0.2494 9.388
## cviolation:cp600layer 4.6014 0.5114 8.998
##
## Correlation of Fixed Effects:
## (Intr) cviltn cp600l
## cviolation 0.819
## cp600layer 0.921 0.891
## cvltn:cp600 0.880 0.932 0.925
## [1] "remove . ~ . - cviolation:cp600layer "## refitting model(s) with ML (instead of REML)## Data: tenseinput
## Models:
## model2: abssum ~ cviolation + cp600layer + (1 + cviolation * cp600layer |
## model2: sub)
## model: abssum ~ cviolation * cp600layer + (1 + cviolation * cp600layer |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 14 3688.9 3760.1 -1830.4 3660.9
## model 15 3667.9 3744.2 -1818.9 3637.9 22.996 1 1.623e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation:cp600layer *** "
## [1] "remove . ~ . - cp600layer "## refitting model(s) with ML (instead of REML)## Data: tenseinput
## Models:
## model2: abssum ~ cviolation + (1 + cviolation * cp600layer | sub)
## model: abssum ~ cviolation + cp600layer + (1 + cviolation * cp600layer |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 13 3687.7 3753.9 -1830.8 3661.7
## model 14 3688.9 3760.1 -1830.4 3660.9 0.845 1 0.358
## [1] "########## Above comparison for cp600layer "
## [1] "remove . ~ . - cviolation "## refitting model(s) with ML (instead of REML)## Data: tenseinput
## Models:
## model2: abssum ~ (1 + cviolation * cp600layer | sub)
## model: abssum ~ cviolation + (1 + cviolation * cp600layer | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 12 3691.2 3752.2 -1833.6 3667.2
## model 13 3687.7 3753.9 -1830.8 3661.7 5.4567 1 0.01949 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation *** "
## $`layer = word`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Control - Violation == 0 -0.83035 0.09871 -8.412 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
##
## $`layer = hidden`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Control - Violation == 0 -5.4317 0.5181 -10.48 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method) # figure for paper
# mapping = aes_string(x="depth", y=dependMeasure, colour = "Condition",linetype="Condition")
drawERP(tenseinput,dependMeasure,"Condition", timeline.y=8,span=.8)
if (!noslopes){ # when doing final version, change noslopes to FALSE
myggsave("img/tense.png",width=figwidth,height=figheight)
}
n4P6cor = computeN4P6correlation(tenseinput)
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.654198792622523"n4P6cor = computeN4P6correlation2(tenseinput,"Control","Violation")
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.339212224579107"CUT AND PASTE INTO LATEX DOC click in box and use right arrow to see whole text
This is an example of the stimuli used in the model.
\begin{tabular}{l|l}
\hline
Condition & Example\\
\hline
Control & a father will drink the beer -s . .\\
\hline
Violation & a father will drink -ed the beer -s . .\\
\hline
\end{tabular}
The maximal model for the data had random slopes + (1 + cviolation*cp600layer |sub)
Mixed Model Results
a main effect of cviolation ($\beta$=3.1, SE=0.27, $\chi^2$(1)=5.46, p=0.019)
no main effect of cp600layer ($\beta$=2.3, SE=0.25, $\chi^2$(1)=0.84, p=0.358)
@@an interaction of cviolation:cp600layer ($\beta$=4.6, SE=0.51, $\chi^2$(1)=23, p<0.001)
Posthoc tests for each layer
@@Posthoc: There was a significant difference for layer = word, diff = 0.8303, t(0)=8.41, p<0.001
@@Posthoc: There was a significant difference for layer = hidden, diff = 5.4317, t(0)=10.48, p<0.001
ratio layer = hidden 6.5418523425268 rev 0.152861903271536WORDCAT (Wassenaar and Hagoort, 2005)
write("wordcat", stderr())
cat.df = subset(adultdata.df, str_detect(Condition,"(CATCONT|CATVIOL)") )
cat.df$Condition=factor(cat.df$Condition,labels=c("Control","Violation"))
verbpos = which(cat.df$wordcat == "VERBI")
wordcatsamelen = which(cat.df$wordcat!=cat.df$wordcat[1])[1]-1
verb.df = cat.df[verbpos,] # verb
postverb.df = cat.df[verbpos + wordcatsamelen,] # position after verb
postverb.df=postverb.df[!is.na(postverb.df$Condition),]
print(head(postverb.df[postverb.df$measure=="target",c("word","wordcat","Condition","Example","tick")])) # check## word wordcat Condition Example
## 113218 -s PL Control the grandma was take -ing the nap -s . .
## 113348 -ed ED Violation the grandma was take -ing the nap -ed . .
## 147265 -s PL Control the girl take -ed the nap -s . .
## 147382 -ed ED Violation the girl take -ed the nap -ed . .
## 160343 -s PL Control a teacher was take -ing the nap -s . .
## 160473 -ed ED Violation a teacher was take -ing the nap -ed . .
## tick
## 113218 8
## 113348 8
## 147265 7
## 147382 7
## 160343 8
## 160473 8showWordOutError("X.ed,X.s,per",c("Condition","measure"),postverb.df)
corrHiddenOutputInputDeriv(postverb.df,"Condition")## Condition
## Control Violation
## 8 8
#mapping = aes(x=depth, y=value, colour = Condition,linetype=Condition)
norawpostverbnotargout.df=removeWordsDrawERP(postverb.df,"value", "Condition")## sub epoch word wordcat contfunc
## 1 s0 100000 -s PL F
## 2 s0 100000 -s PL F
## 3 s0 100000 -s PL F
## 4 s0 100000 -s PL F
## 5 s0 100000 -s PL F
## 6 s0 100000 -s PL F
## Example Condition tick layer
## 1 the grandma was take -ing the nap -s . . Control 8 word
## 2 the grandma was take -ing the nap -s . . Control 8 word
## 3 the grandma was take -ing the nap -s . . Control 8 compress
## 4 the grandma was take -ing the nap -s . . Control 8 compress
## 5 the grandma was take -ing the nap -s . . Control 8 hidden
## 6 the grandma was take -ing the nap -s . . Control 8 hidden
## measure depth variable value
## 1 outderiv 1 abssum 6.3874720
## 2 inputderiv 1 abssum 1.6868865
## 3 outderiv 2 abssum 78.5175257
## 4 inputderiv 2 abssum 0.3688881
## 5 outderiv 3 abssum 9.9624984
## 6 inputderiv 3 abssum 0.8070717interactionText="NO TEXT"
modellist2 = list()
# get example for table
catwordexample = subset(norawpostverbnotargout.df,measure == mainMeasure & layer == "word")
exampleTable = catwordexample[1:4,c("Condition","Example")]
tabletext = kable(exampleTable,format="latex",row.names=FALSE)
print(exampleTable)## Condition Example
## 113220 Control the grandma was take -ing the nap -s . .
## 113350 Violation the grandma was take -ing the nap -ed . .
## 147267 Control the girl take -ed the nap -s . .
## 147384 Violation the girl take -ed the nap -ed . .catinput = subset(norawpostverbnotargout.df,measure == mainMeasure & layer %in% c("word",p600layer) )
if (nlevels(adultdata.df$sub) > 4){ # do mixed with more than 9 subj
catinput$cviolation = ifelse(catinput$Condition=="Violation",0.5,-0.5) # one against two
catinput$cp600layer = ifelse(catinput$layer==p600layer,0.5,-0.5)
# this is the main mixed model with centered variables
randomeff = paste(" + (1 + Condition * layer |sub)")
if (noslopes){ randomeff = " + (1 | sub)" }
formu = as.formula(paste(dependMeasure,"~ cviolation*cp600layer",randomeff))
omnimodel = lmer(formu, catinput,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodel))
modellist2 = modelComparison(omnimodel) # create anova table
# this model is for posthocs
formu2 = as.formula(paste(dependMeasure,"~ Condition*layer",randomeff))
omnimodelFactor = lmer(formu2, catinput)
#print(summary(omnimodelFactor))
model.lsmobj <- lsmeans(omnimodelFactor, ~ Condition | layer)
posthocs = summary(as.glht(pairs(model.lsmobj)))
print(posthocs)
interactionText = printInteraction(posthocs)
difflex = difflex + abs(posthocs$`layer = word`$test$coefficients[[1]])
}## Linear mixed model fit by REML ['lmerMod']
## Formula: abssum ~ cviolation * cp600layer + (1 + Condition * layer | sub)
## Data: catinput
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: 3708.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2596 -0.2688 0.0028 0.1206 7.8185
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.0000 0.0000
## ConditionViolation 0.4319 0.6572 NaN
## layerhidden 0.5353 0.7316 NaN -0.28
## ConditionViolation:layerhidden 4.7542 2.1804 NaN 0.26
## Residual 1.1874 1.0897
##
##
##
##
## -0.55
##
## Number of obs: 1200, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.0250 0.1885 10.740
## cviolation 0.9010 0.4503 2.001
## cp600layer 0.7594 0.2978 2.550
## cviolation:cp600layer 1.7187 0.7009 2.452
##
## Correlation of Fixed Effects:
## (Intr) cviltn cp600l
## cviolation 0.835
## cp600layer 0.800 0.598
## cvltn:cp600 0.707 0.870 0.719
## convergence code: 0
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
##
## [1] "remove . ~ . - cviolation:cp600layer "## refitting model(s) with ML (instead of REML)## Data: catinput
## Models:
## model2: abssum ~ cviolation + cp600layer + (1 + Condition * layer | sub)
## model: abssum ~ cviolation * cp600layer + (1 + Condition * layer | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 14 3728.5 3799.7 -1850.2 3700.5
## model 15 3732.7 3809.1 -1851.4 3702.7 0 1 1
## [1] "########## Above comparison for cviolation:cp600layer "
## [1] "remove . ~ . - cp600layer "## refitting model(s) with ML (instead of REML)## Data: catinput
## Models:
## model2: abssum ~ cviolation + (1 + Condition * layer | sub)
## model: abssum ~ cviolation + cp600layer + (1 + Condition * layer | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 13 3727.5 3793.7 -1850.8 3701.5
## model 14 3728.5 3799.7 -1850.2 3700.5 1.0421 1 0.3073
## [1] "########## Above comparison for cp600layer "
## [1] "remove . ~ . - cviolation "## refitting model(s) with ML (instead of REML)## Data: catinput
## Models:
## model2: abssum ~ (1 + Condition * layer | sub)
## model: abssum ~ cviolation + (1 + Condition * layer | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 12 3732.7 3793.8 -1854.3 3708.7
## model 13 3727.5 3793.7 -1850.8 3701.5 7.1576 1 0.007465 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation *** "
## $`layer = word`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Control - Violation == 0 -0.0417 0.2261 -0.184 0.854
## (Adjusted p values reported -- single-step method)
##
##
## $`layer = hidden`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Control - Violation == 0 -1.7604 0.7747 -2.272 0.0231 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method) # figure for paper
# mapping = aes_string(x="depth", y=dependMeasure, colour = "Condition",linetype="Condition")
drawERP(catinput,dependMeasure,"Condition", timeline.y=4,span=0.8)
if (!noslopes){ # when doing final version, change noslopes to FALSE??n
myggsave("img/cat.png",width=figwidth,height=figheight)
}
n4P6cor = computeN4P6correlation(catinput)
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.400936891846536"n4P6cor = computeN4P6correlation2(catinput,"Control","Violation")
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.451065228107988"CUT AND PASTE INTO LATEX DOC click in box and use right arrow to see whole text
This is an example of the stimuli used in the model.
\begin{tabular}{l|l}
\hline
Condition & Example\\
\hline
Control & the grandma was take -ing the nap -s . .\\
\hline
Violation & the grandma was take -ing the nap -ed . .\\
\hline
Control & the girl take -ed the nap -s . .\\
\hline
Violation & the girl take -ed the nap -ed . .\\
\hline
\end{tabular}
The maximal model for the data had random slopes + (1 + Condition * layer |sub)
Mixed Model Results
a main effect of cviolation ($\beta$=0.9, SE=0.45, $\chi^2$(1)=7.16, p=0.0075)
no main effect of cp600layer ($\beta$=0.76, SE=0.3, $\chi^2$(1)=1.04, p=0.307)
@@no interaction of cviolation:cp600layer ($\beta$=1.7, SE=0.7, $\chi^2$(1)=0, p=1.000)
Posthoc tests for each layer
@@Posthoc: There was no difference for layer = word, p=0.8536
@@Posthoc: There was a significant difference for layer = hidden, diff = 1.7604, t(0)=2.27, p=0.0231
ratio layer = hidden 42.2158273381295 rev 0.0236877982276755SUBCAT Osterhout and Holcomb (1992)
write("subcat", stderr())
sub.df = subset(adultdata.df, str_detect(Condition,"(SUBCONT|SUBVIOL)") )
sub.df$Condition=factor(sub.df$Condition,labels=c("Control","Violation"))
preppos = which(sub.df$word %in% c("near","by"))
prep.df = sub.df[preppos,] # verb
print(head(prep.df[prep.df$measure=="target",c("word","wordcat","Condition","Example","tick")])) # check## word wordcat Condition Example
## 109162 near PREP Control a sister is nap -ing near the boy . .
## 109292 near PREP Violation a sister is take -ing near the boy . .
## 111021 by BY Control the uncle -s walk by tea -s . .
## 111138 by BY Violation the uncle -s drink by tea -s . .
## 112555 by BY Control a grandma jump -ed by the brother -s . .
## 112685 by BY Violation a grandma sip -ed by the brother -s . .
## tick
## 109162 6
## 109292 6
## 111021 5
## 111138 5
## 112555 5
## 112685 5showWordOutError("near,by,the,per",c("Condition","measure"),prep.df)
corrHiddenOutputInputDeriv(prep.df,"Condition")## Condition
## Control Violation
## 8 8
#mapping = aes(x=depth, y=value, colour = Condition,linetype=Condition)
norawprepnotargout.df=removeWordsDrawERP(prep.df,"value","Condition")## sub epoch word wordcat contfunc Example
## 1 s0 100000 near PREP F a sister is nap -ing near the boy . .
## 2 s0 100000 near PREP F a sister is nap -ing near the boy . .
## 3 s0 100000 near PREP F a sister is nap -ing near the boy . .
## 4 s0 100000 near PREP F a sister is nap -ing near the boy . .
## 5 s0 100000 near PREP F a sister is nap -ing near the boy . .
## 6 s0 100000 near PREP F a sister is nap -ing near the boy . .
## Condition tick layer measure depth variable value
## 1 Control 6 word outderiv 1 abssum 6.5775930
## 2 Control 6 word inputderiv 1 abssum 1.6959371
## 3 Control 6 compress outderiv 2 abssum 24.7851990
## 4 Control 6 compress inputderiv 2 abssum 0.2628510
## 5 Control 6 hidden outderiv 3 abssum 5.0414280
## 6 Control 6 hidden inputderiv 3 abssum 0.3902036interactionText="NO TEXT"
modellist2 = list()
# get example for table
catwordexample = subset(norawprepnotargout.df,measure == mainMeasure & layer == "word")
exampleTable = catwordexample[1:4,c("Condition","Example")]
tabletext = kable(exampleTable,format="latex",row.names=FALSE)
print(exampleTable)## Condition Example
## 109164 Control a sister is nap -ing near the boy . .
## 109294 Violation a sister is take -ing near the boy . .
## 111023 Control the uncle -s walk by tea -s . .
## 111140 Violation the uncle -s drink by tea -s . .subinput = subset(norawprepnotargout.df,measure == mainMeasure & layer %in% c("word",p600layer) )
if (nlevels(adultdata.df$sub) > 4){ # do mixed with more than 9 subj
subinput$cviolation = ifelse(subinput$Condition=="Violation",0.5,-0.5) # one against two
subinput$cp600layer = ifelse(subinput$layer==p600layer,0.5,-0.5)
# this is the main mixed model with centered variables
randomeff = paste(" + (1 + cviolation*cp600layer |sub)")
if (noslopes){ randomeff = " + (1 | sub)"}
formu = as.formula(paste(dependMeasure,"~ cviolation*cp600layer",randomeff))
omnimodel = lmer(formu, subinput,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodel))
modellist2 = modelComparison(omnimodel) # create anova table
# this model is for posthocs
formu2 = as.formula(paste(dependMeasure,"~ Condition*layer",randomeff))
omnimodelFactor = lmer(formu2, subinput)
#print(summary(omnimodelFactor))
model.lsmobj <- lsmeans(omnimodelFactor, ~ Condition | layer)
posthocs = summary(as.glht(pairs(model.lsmobj)))
print(posthocs)
interactionText = printInteraction(posthocs)
difflex = difflex + abs(posthocs$`layer = word`$test$coefficients[[1]])
}## Linear mixed model fit by REML ['lmerMod']
## Formula:
## abssum ~ cviolation * cp600layer + (1 + cviolation * cp600layer |
## sub)
## Data: subinput
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: 4511.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.0550 -0.2155 0.0002 0.0798 7.8464
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.5597 0.7481
## cviolation 2.0680 1.4381 0.99
## cp600layer 2.2341 1.4947 1.00 0.99
## cviolation:cp600layer 8.2793 2.8774 0.99 1.00 0.99
## Residual 2.4064 1.5513
## Number of obs: 1200, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.2931 0.2408 13.676
## cviolation 3.6788 0.4635 7.937
## cp600layer 2.9511 0.4811 6.134
## cviolation:cp600layer 6.6299 0.9274 7.149
##
## Correlation of Fixed Effects:
## (Intr) cviltn cp600l
## cviolation 0.958
## cp600layer 0.965 0.958
## cvltn:cp600 0.958 0.963 0.959
## [1] "remove . ~ . - cviolation:cp600layer "## refitting model(s) with ML (instead of REML)## Data: subinput
## Models:
## model2: abssum ~ cviolation + cp600layer + (1 + cviolation * cp600layer |
## model2: sub)
## model: abssum ~ cviolation * cp600layer + (1 + cviolation * cp600layer |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 14 4551.7 4623.0 -2261.8 4523.7
## model 15 4534.7 4611.1 -2252.4 4504.7 18.976 1 1.324e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation:cp600layer *** "
## [1] "remove . ~ . - cp600layer "## refitting model(s) with ML (instead of REML)## Data: subinput
## Models:
## model2: abssum ~ cviolation + (1 + cviolation * cp600layer | sub)
## model: abssum ~ cviolation + cp600layer + (1 + cviolation * cp600layer |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 13 4550.8 4616.9 -2262.4 4524.8
## model 14 4551.7 4623.0 -2261.8 4523.7 1.0741 1 0.3
## [1] "########## Above comparison for cp600layer "
## [1] "remove . ~ . - cviolation "## refitting model(s) with ML (instead of REML)## Data: subinput
## Models:
## model2: abssum ~ (1 + cviolation * cp600layer | sub)
## model: abssum ~ cviolation + (1 + cviolation * cp600layer | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 12 4553.0 4614.1 -2264.5 4529.0
## model 13 4550.8 4616.9 -2262.4 4524.8 4.225 1 0.03983 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation *** "
## $`layer = word`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Control - Violation == 0 -0.3639 0.1267 -2.872 0.00408 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
##
## $`layer = hidden`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Control - Violation == 0 -6.9938 0.9185 -7.614 2.64e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method) # figure for paper
#mapping = aes_string(x="depth", y=dependMeasure, colour = "Condition",linetype="Condition")
drawERP(subinput,dependMeasure,"Condition", timeline.y=4.7)
if (!noslopes){ # when doing final version, change noslopes to FALSE
myggsave("img/subjcat.png",width=figwidth,height=figheight)
}
n4P6cor = computeN4P6correlation(subinput)
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.667942823337718"n4P6cor = computeN4P6correlation2(subinput,"Control","Violation")
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.12828379198526"CUT AND PASTE INTO LATEX DOC click in box and use right arrow to see whole text
This is an example of the stimuli used in the model.
\begin{tabular}{l|l}
\hline
Condition & Example\\
\hline
Control & a sister is nap -ing near the boy . .\\
\hline
Violation & a sister is take -ing near the boy . .\\
\hline
Control & the uncle -s walk by tea -s . .\\
\hline
Violation & the uncle -s drink by tea -s . .\\
\hline
\end{tabular}
The maximal model for the data had random slopes + (1 + cviolation*cp600layer |sub)
Mixed Model Results
a main effect of cviolation ($\beta$=3.7, SE=0.46, $\chi^2$(1)=4.23, p=0.040)
no main effect of cp600layer ($\beta$=3, SE=0.48, $\chi^2$(1)=1.07, p=0.300)
@@an interaction of cviolation:cp600layer ($\beta$=6.6, SE=0.93, $\chi^2$(1)=18.98, p<0.001)
Posthoc tests for each layer
@@Posthoc: There was a significant difference for layer = word, diff = 0.3639, t(0)=2.87, p=0.0041
@@Posthoc: There was a significant difference for layer = hidden, diff = 6.9938, t(0)=7.61, p<0.001
ratio layer = hidden 19.2190162132454 rev 0.052031799593926GARDEN PATH
write("garden", stderr())
garden.df = subset(adultdata.df, str_detect(Condition,"(GARDSC|GARDAMB)") & tick > 4 )
garden.df$Condition=factor(garden.df$Condition,labels=c("Unambiguous","Ambiguous"),levels=c("GARDSC","GARDAMB"))
garden.df$tarword = str_match(garden.df$Example,"(believe|know)(.*?) (the|a) ([A-z]+)( -s)* ([A-z.]+)*")[,7]
garden.df$strlen = sapply(gregexpr("[^ ]+", garden.df$Example), function(x) sum(x > 0))
gardenpos = which(garden.df$word == garden.df$tarword)
gardentar.df = garden.df[gardenpos,] # verb
gardentar.df= subset(gardentar.df,strlen != tick)
gardentar.df$strlen = NULL
gardentar.df$tarword = NULL
print(head(gardentar.df[gardentar.df$measure=="target",c("word","wordcat","Condition","Example","tick")])) # check## word wordcat Condition
## 129988 run VERBI Unambiguous
## 130118 run VERBI Ambiguous
## 145380 nap VERBI Unambiguous
## 145510 nap VERBI Ambiguous
## 147876 walk VERBI Unambiguous
## 148006 walk VERBI Ambiguous
## Example tick
## 129988 the driver believe -ed that the uncle run -ss . . 8
## 130118 the driver believe -ed the uncle run -ss . . 7
## 145380 man -s will know that the girl nap -ss . . 8
## 145510 man -s will know the girl nap -ss . . 7
## 147876 a grandpa know -ed that the sister walk -ss . . 8
## 148006 a grandpa know -ed the sister walk -ss . . 7showWordOutError("is,are,will,was,were,jump,run,walk,nap,per",c("Condition","measure"),gardentar.df)
corrHiddenOutputInputDeriv(gardentar.df,"Condition")## Condition
## Unambiguous Ambiguous
## 8 8
#mapping = aes(x=depth, y=value, colour = Condition,linetype=Condition)
norawgardennotargout.df=removeWordsDrawERP(gardentar.df,"value","Condition")## sub epoch word wordcat contfunc
## 1 s0 100000 run VERBI C
## 2 s0 100000 run VERBI C
## 3 s0 100000 run VERBI C
## 4 s0 100000 run VERBI C
## 5 s0 100000 run VERBI C
## 6 s0 100000 run VERBI C
## Example Condition tick
## 1 the driver believe -ed that the uncle run -ss . . Unambiguous 8
## 2 the driver believe -ed that the uncle run -ss . . Unambiguous 8
## 3 the driver believe -ed that the uncle run -ss . . Unambiguous 8
## 4 the driver believe -ed that the uncle run -ss . . Unambiguous 8
## 5 the driver believe -ed that the uncle run -ss . . Unambiguous 8
## 6 the driver believe -ed that the uncle run -ss . . Unambiguous 8
## layer measure depth variable value
## 1 word outderiv 1 abssum 7.686687
## 2 word inputderiv 1 abssum 1.739810
## 3 compress outderiv 2 abssum 38.730070
## 4 compress inputderiv 2 abssum 0.831418
## 5 hidden outderiv 3 abssum 19.639438
## 6 hidden inputderiv 3 abssum 1.330123interactionText="NO TEXT"
modellist2 = list()
# get example for table
gardenexample = subset(norawgardennotargout.df,measure == mainMeasure & layer == "word")
exampleTable = gardenexample[1:6,c("Condition","Example")]
tabletext = kable(exampleTable,format="latex",row.names=FALSE)
print(exampleTable)## Condition Example
## 129990 Unambiguous the driver believe -ed that the uncle run -ss . .
## 130120 Ambiguous the driver believe -ed the uncle run -ss . .
## 145382 Unambiguous man -s will know that the girl nap -ss . .
## 145512 Ambiguous man -s will know the girl nap -ss . .
## 147878 Unambiguous a grandpa know -ed that the sister walk -ss . .
## 148008 Ambiguous a grandpa know -ed the sister walk -ss . .gardeninput = subset(norawgardennotargout.df,measure == mainMeasure & layer %in% c("word",p600layer) )
if (nlevels(adultdata.df$sub) > 4){ # do mixed with more than 9 subj
gardeninput$cviolation = ifelse(gardeninput$Condition=="Ambiguous",0.5,-0.5) # one against two
gardeninput$cp600layer = ifelse(gardeninput$layer==p600layer,0.5,-0.5)
# this is the main mixed model with centered variables
randomeff = paste(" + (1+ cviolation*cp600layer |sub)") # Condition * layer
if (noslopes){ randomeff = " + (1 | sub)"}
print(paste("************ RandomEffect ->",randomeff))
formu = as.formula(paste(dependMeasure,"~ cviolation*cp600layer",randomeff))
omnimodel = lmer(formu, gardeninput,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodel))
modellist2 = modelComparison(omnimodel) # create anova table
# this model is for posthocs
formu2 = as.formula(paste(dependMeasure,"~ Condition*layer",randomeff))
omnimodelFactor = lmer(formu2, gardeninput)
#print(summary(omnimodelFactor))
model.lsmobj <- lsmeans(omnimodelFactor, ~ Condition | layer)
posthocs = summary(as.glht(pairs(model.lsmobj)))
print(posthocs)
interactionText = printInteraction(posthocs)
difflex = difflex + abs(posthocs$`layer = word`$test$coefficients[[1]])
}## [1] "************ RandomEffect -> + (1+ cviolation*cp600layer |sub)"
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## abssum ~ cviolation * cp600layer + (1 + cviolation * cp600layer |
## sub)
## Data: gardeninput
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: 1698.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.1074 -0.2501 0.0139 0.1583 7.7765
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.1175 0.3428
## cviolation 0.2086 0.4567 0.71
## cp600layer 0.3961 0.6294 1.00 0.75
## cviolation:cp600layer 0.7992 0.8940 0.78 0.99 0.82
## Residual 0.2238 0.4731
## Number of obs: 1200, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.8242 0.1093 16.694
## cviolation 0.7327 0.1470 4.985
## cp600layer -0.1146 0.2009 -0.571
## cviolation:cp600layer 1.1662 0.2879 4.050
##
## Correlation of Fixed Effects:
## (Intr) cviltn cp600l
## cviolation 0.690
## cp600layer 0.981 0.733
## cvltn:cp600 0.757 0.960 0.795
## [1] "remove . ~ . - cviolation:cp600layer "## refitting model(s) with ML (instead of REML)## Data: gardeninput
## Models:
## model2: abssum ~ cviolation + cp600layer + (1 + cviolation * cp600layer |
## model2: sub)
## model: abssum ~ cviolation * cp600layer + (1 + cviolation * cp600layer |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 14 1723.1 1794.3 -847.54 1695.1
## model 15 1714.7 1791.1 -842.36 1684.7 10.365 1 0.001284 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation:cp600layer *** "
## [1] "remove . ~ . - cp600layer "## refitting model(s) with ML (instead of REML)## Data: gardeninput
## Models:
## model2: abssum ~ cviolation + (1 + cviolation * cp600layer | sub)
## model: abssum ~ cviolation + cp600layer + (1 + cviolation * cp600layer |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 13 1730.4 1796.5 -852.18 1704.4
## model 14 1723.1 1794.3 -847.54 1695.1 9.2713 1 0.002328 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cp600layer *** "
## [1] "remove . ~ . - cviolation "## refitting model(s) with ML (instead of REML)## Data: gardeninput
## Models:
## model2: abssum ~ (1 + cviolation * cp600layer | sub)
## model: abssum ~ cviolation + (1 + cviolation * cp600layer | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 12 1729.6 1790.7 -852.79 1705.6
## model 13 1730.4 1796.5 -852.18 1704.4 1.2265 1 0.2681
## [1] "########## Above comparison for cviolation "
## $`layer = word`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Unambiguous - Ambiguous == 0 -0.14965 0.04146 -3.609 0.000307 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
##
## $`layer = hidden`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Unambiguous - Ambiguous == 0 -1.316 0.288 -4.569 4.9e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method) # figure for paper
# mapping = aes_string(x="depth", y=dependMeasure, colour = "Condition",linetype="Condition")
drawERP(gardeninput,dependMeasure, "Condition", timeline.y=2.2)
if (!noslopes){ # when doing final version, change noslopes to FALSE
myggsave("img/garden.png",width=figwidth,height=figheight)
}
n4P6cor = computeN4P6correlation(gardeninput)
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.515584711933041"n4P6cor = computeN4P6correlation2(gardeninput,"Unambiguous","Ambiguous")
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.138453785274273"CUT AND PASTE INTO LATEX DOC click in box and use right arrow to see whole text
UPDATE THIS
This is an example of the stimuli used in the model.
\begin{tabular}{l|l}
\hline
Condition & Example\\
\hline
Unambiguous & the driver believe -ed that the uncle run -ss . .\\
\hline
Ambiguous & the driver believe -ed the uncle run -ss . .\\
\hline
Unambiguous & man -s will know that the girl nap -ss . .\\
\hline
Ambiguous & man -s will know the girl nap -ss . .\\
\hline
Unambiguous & a grandpa know -ed that the sister walk -ss . .\\
\hline
Ambiguous & a grandpa know -ed the sister walk -ss . .\\
\hline
\end{tabular}
The maximal model for the data had random slopes + (1+ cviolation*cp600layer |sub)
Mixed Model Results
no main effect of cviolation ($\beta$=0.73, SE=0.15, $\chi^2$(1)=1.23, p=0.268)
a main effect of cp600layer ($\beta$=-0.11, SE=0.2, $\chi^2$(1)=9.27, p=0.0023)
@@an interaction of cviolation:cp600layer ($\beta$=1.2, SE=0.29, $\chi^2$(1)=10.36, p=0.0013)
Posthoc tests for each layer
@@Posthoc: There was a significant difference for layer = word, diff = 0.1496, t(0)=3.61, p<0.001
@@Posthoc: There was a significant difference for layer = hidden, diff = 1.3158, t(0)=4.57, p<0.001
ratio layer = hidden 8.79545454545454 rev 0.113695090439276Semantic P6 (Kim and Osterhout)
write("semP6", stderr())
semP6.df = subset(adultdata.df, str_detect(Condition,"(PASSCONT|ACTCONT|SEMATTR|NOATTR)") )
semP6.df$Condition = factor(semP6.df$Condition, labels=c("Active Control","Passive Control", "Role Reversal"),levels = c("ACTCONT","PASSCONT","SEMATTR"))
#"NOATTR", "No Attractor",
verbpos = which(semP6.df$wordcat == "VERBT")
wordcatsamelen = which(semP6.df$wordcat!=semP6.df$wordcat[1])[1]-1
verb.df = semP6.df[verbpos,] # verb
postverb.df = semP6.df[verbpos + wordcatsamelen,] # position after verb
postverb.df=postverb.df[!is.na(postverb.df$Condition),]
print(head(postverb.df[postverb.df$measure=="target",c("word","wordcat","Condition","Example","tick")])) # check## word wordcat Condition
## 109409 -par PAR Passive Control
## 109539 -ing ING Active Control
## 109656 -ing ING Role Reversal
## 111255 -par PAR Passive Control
## 111385 -ing ING Active Control
## 111502 -ing ING Role Reversal
## Example tick
## 109409 the tea is sip -par by a friend . . 5
## 109539 the friend is sip -ing a tea . . 5
## 109656 the tea is sip -ing a friend . . 5
## 111255 a wine is drink -par by the driver . . 5
## 111385 a driver is drink -ing the wine . . 5
## 111502 a wine is drink -ing the driver . . 5showWordOutError("the,a,X.ss,X.ing,X.ed,X.par",c("Condition","measure"),postverb.df)
showWordOutError("the,a,X.ss,X.ing,X.ed,X.par",c("Condition","measure"),verb.df)
computeClozeERPCorr(postverb.df,dependMeasure,"word")## [1] "Correlation output prediction and abssum error of word is -0.999999999999951"computeClozeERPCorr(postverb.df,dependMeasure,"hidden")## [1] "Correlation output prediction and abssum error of hidden is -0.895965039051418"corrHiddenOutputInputDeriv(postverb.df,"Condition")## Condition
## Active Control Passive Control Role Reversal
## 8 8 8
#mapping = aes_string(x="depth", y="value", colour = "Condition",linetype="Condition")
norawpostverbnotargout.df=removeWordsDrawERP(postverb.df,"value","Condition")## sub epoch word wordcat contfunc Example
## 1 s0 100000 -par PAR F the tea is sip -par by a friend . .
## 2 s0 100000 -par PAR F the tea is sip -par by a friend . .
## 3 s0 100000 -par PAR F the tea is sip -par by a friend . .
## 4 s0 100000 -par PAR F the tea is sip -par by a friend . .
## 5 s0 100000 -par PAR F the tea is sip -par by a friend . .
## 6 s0 100000 -par PAR F the tea is sip -par by a friend . .
## Condition tick layer measure depth variable value
## 1 Passive Control 5 word outderiv 1 abssum 1.0011830
## 2 Passive Control 5 word inputderiv 1 abssum 0.0023622
## 3 Passive Control 5 compress outderiv 2 abssum 0.0486240
## 4 Passive Control 5 compress inputderiv 2 abssum 0.0024007
## 5 Passive Control 5 hidden outderiv 3 abssum 0.0407550
## 6 Passive Control 5 hidden inputderiv 3 abssum 0.0035750norawverbnotargout.df=removeWordsDrawERP(verb.df,"value","Condition")## sub epoch word wordcat contfunc Example
## 1 s0 100000 sip VERBT C the tea is sip -par by a friend . .
## 2 s0 100000 sip VERBT C the tea is sip -par by a friend . .
## 3 s0 100000 sip VERBT C the tea is sip -par by a friend . .
## 4 s0 100000 sip VERBT C the tea is sip -par by a friend . .
## 5 s0 100000 sip VERBT C the tea is sip -par by a friend . .
## 6 s0 100000 sip VERBT C the tea is sip -par by a friend . .
## Condition tick layer measure depth variable value
## 1 Passive Control 4 word outderiv 1 abssum 40.5763800
## 2 Passive Control 4 word inputderiv 1 abssum 1.9507104
## 3 Passive Control 4 compress outderiv 2 abssum 31.8691877
## 4 Passive Control 4 compress inputderiv 2 abssum 0.9560113
## 5 Passive Control 4 hidden outderiv 3 abssum 25.9320234
## 6 Passive Control 4 hidden inputderiv 3 abssum 1.8636249interactionText="NO TEXT"
modellist2 = list()
#norawpostverbnotargout.df=subset(norawpostverbnotargout.df,Condition != "No Attractor")
#norawpostverbnotargout.df=norawpostverbnotargout.df
# get example for table
semP6exampleword = subset(norawpostverbnotargout.df,measure == mainMeasure & layer == "word")
exampleTable = semP6exampleword[1:12,c("Condition","Example")]
tabletext = kable(exampleTable,format="latex",row.names=FALSE)
print(exampleTable)## Condition Example
## 109411 Passive Control the tea is sip -par by a friend . .
## 109541 Active Control the friend is sip -ing a tea . .
## 109658 Role Reversal the tea is sip -ing a friend . .
## 111257 Passive Control a wine is drink -par by the driver . .
## 111387 Active Control a driver is drink -ing the wine . .
## 111504 Role Reversal a wine is drink -ing the driver . .
## 112193 Passive Control the water is drink -par by the clerk . .
## 112323 Active Control the clerk is drink -ing the water . .
## 112440 Role Reversal the water is drink -ing the clerk . .
## 112817 Passive Control a tea is taste -par by the aunt . .
## 112947 Active Control a aunt is taste -ing the tea . .
## 113064 Role Reversal a tea is taste -ing the aunt . .semP6input = subset(norawpostverbnotargout.df,measure == mainMeasure & layer %in% c("word",p600layer) )
semP6inputverb = subset(norawverbnotargout.df,measure == mainMeasure & layer %in% c("word",p600layer) )
if (nlevels(adultdata.df$sub) > 4){ # do mixed with more than 9 subj
options(contrasts=c("contr.treatment", "contr.treatment"))
# semP6input$Condition = factor(semP6input$Condition,levels=c("Active Control", "Passive Control", "No Attractor", "Semantic Attractor"))
semP6input$cpassive = 0
semP6input$cpassive[semP6input$Condition=="Passive Control"]= 0.5
semP6input$cpassive[semP6input$Condition=="Active Control"]= -0.5
semP6input$semattvsgram = 0
semP6input$semattvsgram[semP6input$Condition=="Role Reversal"] = 1
semP6input$semattvsgram[semP6input$Condition%in% c("Passive Control","Active Control")] = -0.5
semP6input$cp600layer = ifelse(semP6input$layer==p600layer,0.5,-0.5)
# this is the main mixed model with centered variables
randomeff = paste(" + (1 + cp600layer |sub)") # Condition * layer
if (noslopes){ randomeff = " + (1 | sub)"}
formu = as.formula(paste(dependMeasure,"~ cp600layer*cpassive + cp600layer*semattvsgram ",randomeff))
omnimodel = lmer(formu, semP6input,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodel))
modellist2 = modelComparison(omnimodel) # create anova table
maincond = anova(modellist2[[3]],modellist2[[5]])
condlayer = anova(modellist2[[1]],modellist2[[3]])
# this model is for posthocs
semP6input$semattvsgramf = factor(semP6input$semattvsgram,labels=c("ActPas","RoleRev"))
formu2 = as.formula(paste(dependMeasure,"~ layer*semattvsgramf",randomeff))
omnimodelFactor = lmer(formu2, semP6input)
#print(summary(omnimodelFactor))
model.lsmobj <- lsmeans(omnimodelFactor, ~ semattvsgramf | layer)
posthocs = summary(as.glht(pairs(model.lsmobj)))
print(posthocs)
interactionText = printInteraction(posthocs)
difflexSemP6 = abs(posthocs$`layer = word`$test$coefficients[[1]])
avedifflex = difflex/5
print(paste("Ratio of lexical sem P6 vs other P6",difflexSemP6/avedifflex))
# semP6inputverb$cpassive = 0
# semP6inputverb$cpassive[semP6inputverb$Condition=="Passive Control"]= 0.5
# semP6inputverb$cpassive[semP6inputverb$Condition=="Active Control"]= -0.5
# semP6inputverb$semattvsgram = 0
# semP6inputverb$semattvsgram[semP6inputverb$Condition=="Role Reversal"] = 1
# semP6inputverb$semattvsgram[semP6inputverb$Condition%in% c("Passive Control","Active Control")] = -0.5
# semP6inputverb$cp600layer = ifelse(semP6inputverb$layer==p600layer,0.5,-0.5)
#
# # this is the main mixed model with centered variables
# randomeff = paste(" + (1 + cp600layer |sub)") # Condition * layer
# if (noslopes){ randomeff = " + (1 | sub)"}
# formu = as.formula(paste(dependMeasure,"~ cp600layer*cpassive + cp600layer*semattvsgram ",randomeff))
# omnimodel2 = lmer(formu, semP6inputverb,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
# print(summary(omnimodel2))
# modellist2 = modelComparison(omnimodel) # create anova table
#
# maincond = anova(modellist2[[3]],modellist2[[5]])
# condlayer = anova(modellist2[[1]],modellist2[[3]])
#
# # this model is for posthocs
# formu2 = as.formula(paste(dependMeasure,"~ layer*Condition",randomeff))
# omnimodelFactor = lmer(formu2, semP6inputverb)
# #print(summary(omnimodelFactor))
# model.lsmobj <- lsmeans(omnimodelFactor, ~ Condition | layer)
# posthocs = summary(as.glht(pairs(model.lsmobj)))
# print(posthocs)
# interactionText = printInteraction(posthocs)
}## Linear mixed model fit by REML ['lmerMod']
## Formula: abssum ~ cp600layer * cpassive + cp600layer * semattvsgram +
## (1 + cp600layer | sub)
## Data: semP6input
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: 5325.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3154 -0.0211 -0.0001 0.0714 5.5308
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.2442 0.4942
## cp600layer 0.9793 0.9896 1.00
## Residual 1.0926 1.0453
## Number of obs: 1800, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.66508 0.15821 10.525
## cp600layer 1.99271 0.31680 6.290
## cpassive 0.07462 0.06035 1.237
## semattvsgram 3.24257 0.03484 93.063
## cp600layer:cpassive 0.10791 0.12070 0.894
## cp600layer:semattvsgram 3.86122 0.06969 55.409
##
## Correlation of Fixed Effects:
## (Intr) cp600l cpassv smttvs cp600lyr:c
## cp600layer 0.976
## cpassive 0.000 0.000
## semattvsgrm 0.000 0.000 0.000
## cp600lyr:cp 0.000 0.000 0.000 0.000
## cp600lyr:sm 0.000 0.000 0.000 0.000 0.000
## [1] "remove . ~ . - cp600layer:semattvsgram "## refitting model(s) with ML (instead of REML)## Data: semP6input
## Models:
## model2: abssum ~ cp600layer + cpassive + semattvsgram + (1 + cp600layer |
## model2: sub) + cp600layer:cpassive
## model: abssum ~ cp600layer * cpassive + cp600layer * semattvsgram +
## model: (1 + cp600layer | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 9 7114.7 7164.1 -3548.3 7096.7
## model 10 5325.6 5380.5 -2652.8 5305.6 1791.1 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cp600layer:semattvsgram *** "
## [1] "remove . ~ . - cp600layer:cpassive "## refitting model(s) with ML (instead of REML)## Data: semP6input
## Models:
## model2: abssum ~ cp600layer + cpassive + semattvsgram + (1 + cp600layer |
## model2: sub)
## model: abssum ~ cp600layer + cpassive + semattvsgram + (1 + cp600layer |
## model: sub) + cp600layer:cpassive
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 8 7113.0 7156.9 -3548.5 7097.0
## model 9 7114.7 7164.1 -3548.3 7096.7 0.2947 1 0.5872
## [1] "########## Above comparison for cp600layer:cpassive "
## [1] "remove . ~ . - semattvsgram "## refitting model(s) with ML (instead of REML)## Data: semP6input
## Models:
## model2: abssum ~ cp600layer + cpassive + (1 + cp600layer | sub)
## model: abssum ~ cp600layer + cpassive + semattvsgram + (1 + cp600layer |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 7 8943.4 8981.9 -4464.7 8929.4
## model 8 7113.0 7156.9 -3548.5 7097.0 1832.4 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for semattvsgram *** "
## [1] "remove . ~ . - cpassive "## refitting model(s) with ML (instead of REML)## Data: semP6input
## Models:
## model2: abssum ~ cp600layer + (1 + cp600layer | sub)
## model: abssum ~ cp600layer + cpassive + (1 + cp600layer | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 6 8941.6 8974.6 -4464.8 8929.6
## model 7 8943.4 8981.9 -4464.7 8929.4 0.2025 1 0.6527
## [1] "########## Above comparison for cpassive "
## [1] "remove . ~ . - cp600layer "## refitting model(s) with ML (instead of REML)## Data: semP6input
## Models:
## model2: abssum ~ (1 + cp600layer | sub)
## model: abssum ~ cp600layer + (1 + cp600layer | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 5 8955.9 8983.4 -4472.9 8945.9
## model 6 8941.6 8974.6 -4464.8 8929.6 16.242 1 5.574e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cp600layer *** "## refitting model(s) with ML (instead of REML)
## refitting model(s) with ML (instead of REML)## $`layer = word`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## ActPas - RoleRev == 0 -1.96794 0.07392 -26.62 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
##
## $`layer = hidden`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## ActPas - RoleRev == 0 -7.75976 0.07392 -105 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
##
## [1] "Ratio of lexical sem P6 vs other P6 5.08163295491971" # figure for paper
semP6input[semP6input$Condition == "Passive Control",dependMeasure]= 0.1+semP6input[semP6input$Condition == "Passive Control",dependMeasure]
# mapping = aes_string(x="depth", y=dependMeasure, colour = "Condition",linetype="Condition")
p = drawERP(semP6input,dependMeasure,"Condition", timeline.y=6)
p = p + scale_linetype_manual(values=c("dotted","dashed","solid")) #
p
if (!noslopes){ # when doing final version, change noslopes to FALSE
myggsave("img/semP6.png",width=figwidth,height=figheight)
}
# mapping = aes_string(x="depth", y=dependMeasure, colour = "Condition",linetype="Condition")
# drawERP(semP6inputverb,mapping, timeline.y=6)
n4P6cor = computeN4P6correlation(semP6input)
print(paste("N400/P600 correlation",n4P6cor))## [1] "N400/P600 correlation -0.894347383940244"CUT AND PASTE INTO LATEX DOC click in box and use right arrow to see whole text
This is an example of the stimuli used in the model.
\begin{tabular}{l|l}
\hline
Condition & Example\\
\hline
Passive Control & the tea is sip -par by a friend . .\\
\hline
Active Control & the friend is sip -ing a tea . .\\
\hline
Role Reversal & the tea is sip -ing a friend . .\\
\hline
Passive Control & a wine is drink -par by the driver . .\\
\hline
Active Control & a driver is drink -ing the wine . .\\
\hline
Role Reversal & a wine is drink -ing the driver . .\\
\hline
Passive Control & the water is drink -par by the clerk . .\\
\hline
Active Control & the clerk is drink -ing the water . .\\
\hline
Role Reversal & the water is drink -ing the clerk . .\\
\hline
Passive Control & a tea is taste -par by the aunt . .\\
\hline
Active Control & a aunt is taste -ing the tea . .\\
\hline
Role Reversal & a tea is taste -ing the aunt . .\\
\hline
\end{tabular}
The maximal model for the data had random slopes + (1 + cp600layer |sub)
main effect $\chi^2$(2)=1832.65, p<0.001
interaction $\chi^2$(2)=1791.41, p<0.001
Mixed Model Results
a main effect of cp600layer ($\beta$=2, SE=0.32, $\chi^2$(1)=16.24, p<0.001)
no main effect of cpassive ($\beta$=0.075, SE=0.06, $\chi^2$(1)=0.2, p=0.653)
a main effect of semattvsgram ($\beta$=3.2, SE=0.035, $\chi^2$(1)=1832.45, p<0.001)
no interaction of cp600layer:cpassive ($\beta$=0.11, SE=0.12, $\chi^2$(1)=0.29, p=0.587)
@@an interaction of cp600layer:semattvsgram ($\beta$=3.9, SE=0.07, $\chi^2$(1)=1791.12, p<0.001)
Posthoc tests for each layer
@@Posthoc: There was a significant difference for layer = word, diff = 1.9679, t(0)=26.62, p<0.001
@@Posthoc: There was a significant difference for layer = hidden, diff = 7.7598, t(0)=104.98, p<0.001
ratio layer = hidden 3.9431881701306 rev 0.253601896956107WORDCAT DEVELOP
write("dev wordcat", stderr())
cat.df = subset(develop.df,str_detect(Condition,"(CATCONT|CATVIOL)") & ! measure %in% c("error", "myoutputder","myinputder","outderiv"))
cat.df$Condition=factor(cat.df$Condition,labels=c("Control","Violation"))
verbpos = which(cat.df$wordcat == "VERBI")
wordcatsamelen = which(cat.df$wordcat!=cat.df$wordcat[1])[1]-1
verb.df = cat.df[verbpos,] # verb
postverb.df = cat.df[verbpos + wordcatsamelen,] # position after verb
postverb.df=postverb.df[!is.na(postverb.df$Condition),]
postverb.df=postverb.df[!is.na(postverb.df$Condition),]
print(head(postverb.df[postverb.df$measure=="target",c("word","wordcat","Condition","Example","tick","epoch")])) # check## word wordcat Condition Example
## 26820 -s PL Control the grandma was take -ing the nap -s . .
## 26950 -ed ED Violation the grandma was take -ing the nap -ed . .
## 60867 -s PL Control the girl take -ed the nap -s . .
## 60984 -ed ED Violation the girl take -ed the nap -ed . .
## 73945 -s PL Control a teacher was take -ing the nap -s . .
## 74075 -ed ED Violation a teacher was take -ing the nap -ed . .
## tick epoch
## 26820 8 10000
## 26950 8 10000
## 60867 7 10000
## 60984 7 10000
## 73945 8 10000
## 74075 8 10000postverb2.df = subset(postverb.df, layer == "word")
postverb2.df$X.s=abs(postverb2.df$X.s)
postverb2.df$X.ss=abs(postverb2.df$X.ss)
postverb2.df$per=abs(postverb2.df$per)
p1=showWordOutError("X.ed,X.s,per",c("Condition","epoch","measure"),postverb2.df)
p1=p1+theme_bw()+ylab("Activation/Error")
#p1 =p1 + scale_fill_discrete(name="Word Unit",labels=c("3rdSingular", "Plural", "End of sentence"))
p1 = p1+ theme(legend.position="bottom")
p1
#myggsave("img/devoutput.png",width=figwidth,height=figheight)corrHiddenOutputInputDeriv(postverb.df,"Condition","epoch")## Condition
## Control Violation
## 8 8
epochnum = which(names(postverb.df)=="epoch")
postverb2.df = postverb.df[,c(epochnum,1:epochnum-1,(epochnum+1):length(postverb.df))]
mapping = aes(x=depth, y=value, colour = Condition,linetype=Condition)
norawpostverbnotargout.df=removeWordsDrawERP(postverb.df,"value", "Condition", "epoch")## sub epoch word wordcat contfunc
## 1 s0 10000 -s PL F
## 2 s0 10000 -s PL F
## 3 s0 10000 -s PL F
## 4 s0 10000 -ed ED F
## 5 s0 10000 -ed ED F
## 6 s0 10000 -ed ED F
## Example Condition tick layer
## 1 the grandma was take -ing the nap -s . . Control 8 word
## 2 the grandma was take -ing the nap -s . . Control 8 compress
## 3 the grandma was take -ing the nap -s . . Control 8 hidden
## 4 the grandma was take -ing the nap -ed . . Violation 8 word
## 5 the grandma was take -ing the nap -ed . . Violation 8 compress
## 6 the grandma was take -ing the nap -ed . . Violation 8 hidden
## measure depth variable value
## 1 inputderiv 1 abssum 0.9018836
## 2 inputderiv 2 abssum 1.2312222
## 3 inputderiv 3 abssum 1.8882032
## 4 inputderiv 1 abssum 1.9999993
## 5 inputderiv 2 abssum 2.5949130
## 6 inputderiv 3 abssum 3.0311034interactionText="NO TEXT"
modellist2 = list()
childtime = 30000
certainepochs.df = subset(norawpostverbnotargout.df, epoch %in% c(childtime,100000))
#certainepochs.df = norawpostverbnotargout.df
# get example for table
catwordexample = subset(certainepochs.df,measure == mainMeasure & layer == "word")
exampleTable = catwordexample[1:4,c("Condition","Example")]
tabletext = kable(exampleTable,format="latex",row.names=FALSE)
print(exampleTable)## Condition Example
## 113220 Control the grandma was take -ing the nap -s . .
## 113350 Violation the grandma was take -ing the nap -ed . .
## 147267 Control the girl take -ed the nap -s . .
## 147384 Violation the girl take -ed the nap -ed . .catinput = subset(certainepochs.df,measure == mainMeasure & layer %in% c("word",p600layer) )
if (nlevels(adultdata.df$sub) > 4){ # do mixed with more than 9 subj
catinput$cviolation = ifelse(catinput$Condition=="Violation",0.5,-0.5) # one against two
catinput$cword = ifelse(catinput$layer=="word",0.5,-0.5)
catinput$cepoch = scale(catinput$epoch)
catinput$sub = factor(catinput$sub)
# this is the main mixed model with centered variables
randomeff = paste(" + (1 + Condition |sub)") # Condition * epoch
if (noslopes){ randomeff = " + (1 | sub)"}
formu = as.formula(paste(dependMeasure,"~ cviolation*cepoch*cword",randomeff))
omnimodel = lmer(formu, catinput,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodel))
modellist2 = modelComparison(omnimodel) # create anova table
catinputword = subset(catinput,layer == "word")
catinputword$cviolation = ifelse(catinputword$Condition=="Violation",0.5,-0.5) # one against two
catinputword$cword = ifelse(catinputword$layer=="word",0.5,-0.5)
catinputword$cepoch = scale(catinputword$epoch)
catinputword$sub = factor(catinputword$sub)
# this is the main mixed model with centered variables
randomeff = paste(" + (1 + Condition |sub)") # Condition * epoch
if (noslopes){ randomeff = " + (1 | sub)"}
formu = as.formula(paste(dependMeasure,"~ cviolation*cepoch*cword",randomeff))
omnimodel = lmer(formu, catinputword,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodel))
modellist2 = modelComparison(omnimodel) # create anova table
# this model is for posthocs
catinputword$depoch = factor(catinputword$epoch)
formu2 = as.formula(paste(dependMeasure,"~ Condition*depoch",randomeff))
omnimodelFactor = lmer(formu2, catinputword)
#print(summary(omnimodelFactor))
model.lsmobj <- lsmeans(omnimodelFactor, ~ Condition | depoch)
posthocs = summary(as.glht(pairs(model.lsmobj)))
print(posthocs)
interactionText = printInteraction(posthocs)
}## Linear mixed model fit by REML ['lmerMod']
## Formula: abssum ~ cviolation * cepoch * cword + (1 + Condition | sub)
## Data: catinput
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: 8791.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9890 -0.5282 -0.0220 0.3025 8.1033
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.2620 0.5118
## ConditionViolation 0.6876 0.8292 -0.60
## Residual 2.2118 1.4872
## Number of obs: 2400, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.41650 0.13654 17.698
## cviolation 1.31298 0.26915 4.878
## cepoch -0.39155 0.03036 -12.895
## cword -1.49342 0.06072 -24.597
## cviolation:cepoch -0.41204 0.06073 -6.785
## cviolation:cword -1.99491 0.12143 -16.428
## cepoch:cword 0.73421 0.06073 12.090
## cviolation:cepoch:cword 0.27632 0.12146 2.275
##
## Correlation of Fixed Effects:
## (Intr) cviltn cepoch cword cvltn:cp cvltn:cw cpch:c
## cviolation 0.237
## cepoch 0.000 0.000
## cword 0.000 0.000 0.000
## cviltn:cpch 0.000 0.000 0.000 0.000
## cviltn:cwrd 0.000 0.000 0.000 0.000 0.000
## cepoch:cwrd 0.000 0.000 0.000 0.000 0.000 0.000
## cvltn:cpch: 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## [1] "remove . ~ . - cviolation:cepoch:cword "## refitting model(s) with ML (instead of REML)## Data: catinput
## Models:
## model2: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub) +
## model2: cviolation:cepoch + cviolation:cword + cepoch:cword
## model: abssum ~ cviolation * cepoch * cword + (1 + Condition | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 11 8794.4 8858.0 -4386.2 8772.4
## model 12 8791.2 8860.6 -4383.6 8767.2 5.1834 1 0.0228 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation:cepoch:cword *** "
## [1] "remove . ~ . - cepoch:cword "## refitting model(s) with ML (instead of REML)## Data: catinput
## Models:
## model2: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub) +
## model2: cviolation:cepoch + cviolation:cword
## model: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub) +
## model: cviolation:cepoch + cviolation:cword + cepoch:cword
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 10 8934.3 8992.1 -4457.1 8914.3
## model 11 8794.4 8858.0 -4386.2 8772.4 141.91 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cepoch:cword *** "
## [1] "remove . ~ . - cviolation:cword "## refitting model(s) with ML (instead of REML)## Data: catinput
## Models:
## model2: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub) +
## model2: cviolation:cepoch
## model: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub) +
## model: cviolation:cepoch + cviolation:cword
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 9 9173.9 9226.0 -4578.0 9155.9
## model 10 8934.3 8992.1 -4457.1 8914.3 241.67 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation:cword *** "
## [1] "remove . ~ . - cviolation:cepoch "## refitting model(s) with ML (instead of REML)## Data: catinput
## Models:
## model2: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub)
## model: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub) +
## model: cviolation:cepoch
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 8 9210.8 9257.1 -4597.4 9194.8
## model 9 9173.9 9226.0 -4578.0 9155.9 38.879 1 4.51e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation:cepoch *** "
## [1] "remove . ~ . - cword "## refitting model(s) with ML (instead of REML)## Data: catinput
## Models:
## model2: abssum ~ cviolation + cepoch + (1 + Condition | sub)
## model: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 7 9668.3 9708.8 -4827.1 9654.3
## model 8 9210.8 9257.1 -4597.4 9194.8 459.45 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cword *** "
## [1] "remove . ~ . - cepoch "## refitting model(s) with ML (instead of REML)## Data: catinput
## Models:
## model2: abssum ~ cviolation + (1 + Condition | sub)
## model: abssum ~ cviolation + cepoch + (1 + Condition | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 6 9778.4 9813.1 -4883.2 9766.4
## model 7 9668.3 9708.8 -4827.1 9654.3 112.15 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cepoch *** "
## [1] "remove . ~ . - cviolation "## refitting model(s) with ML (instead of REML)## Data: catinput
## Models:
## model2: abssum ~ (1 + Condition | sub)
## model: abssum ~ cviolation + (1 + Condition | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 5 9789.4 9818.3 -4889.7 9779.4
## model 6 9778.4 9813.1 -4883.2 9766.4 12.931 1 0.0003231 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation *** "## fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients## Linear mixed model fit by REML ['lmerMod']
## Formula: abssum ~ cviolation * cepoch * cword + (1 + Condition | sub)
## Data: catinputword
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: 672
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5396 -0.5152 0.0633 0.5340 2.8717
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.02381 0.1543
## ConditionViolation 0.17298 0.4159 -0.56
## Residual 0.09528 0.3087
## Number of obs: 1200, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.669794 0.056770 29.413
## cviolation 0.315527 0.132723 2.377
## cepoch -0.024448 0.008915 -2.742
## cviolation:cepoch -0.273938 0.017829 -15.365
##
## Correlation of Fixed Effects:
## (Intr) cviltn cepoch
## cviolation 0.675
## cepoch 0.000 0.000
## cviltn:cpch 0.000 0.000 0.000
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients
## [1] "remove . ~ . - cviolation:cepoch:cword "## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients
## refitting model(s) with ML (instead of REML)## Data: catinputword
## Models:
## model2: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub) +
## model2: cviolation:cepoch + cviolation:cword + cepoch:cword
## model: abssum ~ cviolation * cepoch * cword + (1 + Condition | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 8 667.38 708.1 -325.69 651.38
## model 8 667.38 708.1 -325.69 651.38 0 0 1
## [1] "########## Above comparison for cviolation:cepoch:cword "
## [1] "remove . ~ . - cepoch:cword "## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## refitting model(s) with ML (instead of REML)## Data: catinputword
## Models:
## model2: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub) +
## model2: cviolation:cepoch + cviolation:cword
## model: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub) +
## model: cviolation:cepoch + cviolation:cword + cepoch:cword
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 8 667.38 708.1 -325.69 651.38
## model 8 667.38 708.1 -325.69 651.38 0 0 1
## [1] "########## Above comparison for cepoch:cword "
## [1] "remove . ~ . - cviolation:cword "## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## refitting model(s) with ML (instead of REML)## Data: catinputword
## Models:
## model2: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub) +
## model2: cviolation:cepoch
## model: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub) +
## model: cviolation:cepoch + cviolation:cword
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 8 667.38 708.1 -325.69 651.38
## model 8 667.38 708.1 -325.69 651.38 0 0 1
## [1] "########## Above comparison for cviolation:cword "
## [1] "remove . ~ . - cviolation:cepoch "## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## refitting model(s) with ML (instead of REML)## Data: catinputword
## Models:
## model2: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub)
## model: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub) +
## model: cviolation:cepoch
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 7 880.91 916.54 -433.46 866.91
## model 8 667.38 708.10 -325.69 651.38 215.53 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation:cepoch *** "
## [1] "remove . ~ . - cword "## refitting model(s) with ML (instead of REML)## Data: catinputword
## Models:
## model2: abssum ~ cviolation + cepoch + (1 + Condition | sub)
## model: abssum ~ cviolation + cepoch + cword + (1 + Condition | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 7 880.91 916.54 -433.46 866.91
## model 7 880.91 916.54 -433.46 866.91 0 0 1
## [1] "########## Above comparison for cword "
## [1] "remove . ~ . - cepoch "## refitting model(s) with ML (instead of REML)## Data: catinputword
## Models:
## model2: abssum ~ cviolation + (1 + Condition | sub)
## model: abssum ~ cviolation + cepoch + (1 + Condition | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 6 885.17 915.71 -436.58 873.17
## model 7 880.91 916.54 -433.46 866.91 6.2593 1 0.01235 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cepoch *** "
## [1] "remove . ~ . - cviolation "## refitting model(s) with ML (instead of REML)## Data: catinputword
## Models:
## model2: abssum ~ (1 + Condition | sub)
## model: abssum ~ cviolation + (1 + Condition | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 5 888.04 913.49 -439.02 878.04
## model 6 885.17 915.71 -436.58 873.17 4.8733 1 0.02727 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation *** "
## $`depoch = 30000`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Control - Violation == 0 -0.5894 0.1339 -4.401 1.08e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
##
## $`depoch = 100000`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Control - Violation == 0 -0.0417 0.1339 -0.311 0.755
## (Adjusted p values reported -- single-step method)form = as.formula(paste(dependMeasure," ~ Condition + epoch + layer",sep=""))
meandf = aggregate(form, catinput, mean)
mapping = aes_string(x="epoch", y=dependMeasure, fill = "Condition")
meandf$epoch= factor(meandf$epoch)
p=ggplot(meandf, mapping)
p=p+geom_bar(position="dodge",stat="identity")
p = p + scale_fill_brewer(palette = "Set1")
p=p+facet_wrap(~ layer)
p=p + theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"))
p
catinputword = subset(catinput, layer == "word")
form = as.formula(paste(dependMeasure," ~ Condition + epoch ",sep=""))
meandf = aggregate(form, catinputword, mean)
mapping = aes_string(x="epoch", y=dependMeasure, fill = "Condition")
meandf$epoch= factor(meandf$epoch,labels = c(paste(childtime,"(child)",sep=""),"100000(adult)"))
p=ggplot(meandf, mapping)
p=p+geom_bar(position="dodge",stat="identity")+theme_bw()
#p = p + scale_fill_brewer(palette = "Set1")
p = p + scale_fill_grey(start = 0.25, end = 0.75)
p = p +ylab("Sum Abs. Error")+xlab("Testing Epoch")
p=p + theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"))
p
if (!noslopes){ # when doing final version, change noslopes to FALSE??n
myggsave("img/catdevelop.png",width=figwidth,height=figheight)
}CUT AND PASTE INTO LATEX DOC click in box and use right arrow to see whole text
This is an example of the stimuli used in the model.
Coulson et al
# This is Coulson et al
write("Coulson et al", stderr())
coulagree.df = subset(couldf, str_detect(Condition,"(SING|PLUR)") )
verbpos = which(coulagree.df$wordcat == "VERBT")
wordcatsamelen = which(coulagree.df$wordcat!=coulagree.df$wordcat[1])[1]-1
verb.df = coulagree.df[verbpos,] # verb
postverbcoul.df = coulagree.df[verbpos + wordcatsamelen,] # position after verb
postverbcoul.df=postverbcoul.df[!is.na(postverbcoul.df$Condition),]
postverbcoul.df = cbind(Probability="",Number="",Agreement="",postverbcoul.df)
cond = str_split_fixed(postverbcoul.df$Condition,";",2)
postverbcoul.df$Number = cond[,1]
postverbcoul.df$Agreement = cond[,2]
postverbcoul.df$Number = factor(postverbcoul.df$Number,labels=c("Singular"))
postverbcoul.df$Agreement = factor(postverbcoul.df$Agreement,labels=c("Control","Violation"))
postverbcoul.df$Probability=NA
postverbcoul.df$Probability[postverbcoul.df$Agreement=="Violation" & str_detect(postverbcoul.df$coul,"vcoul1")] = "Improbable"
postverbcoul.df$Probability[postverbcoul.df$Agreement=="Violation" & str_detect(postverbcoul.df$coul,"vcoul2")] = "Probable"
postverbcoul.df$Probability[postverbcoul.df$Agreement=="Control" & str_detect(postverbcoul.df$coul,"vcoul2")] = "Improbable"
postverbcoul.df$Probability[postverbcoul.df$Agreement=="Control" & str_detect(postverbcoul.df$coul,"vcoul1")] = "Probable"
postverbcoul.df$Probability=factor(postverbcoul.df$Probability)
xtabs(~ Agreement +Probability,postverbcoul.df)## Probability
## Agreement Improbable Probable
## Control 600 2400
## Violation 600 2400print(head(postverbcoul.df[postverbcoul.df$measure=="target",c("word","wordcat","Condition","Example","tick")])) # check## word wordcat Condition Example tick
## 157 -ss SS SING;CONT the boy take -ss a stick . . 4
## 391 wine NOUNI PLUR;CONT the father -s taste wine -s . . 5
## 612 a DET SING;VIOL the mother push a kite . . 4
## 846 the DET PLUR;CONT the wife -s take the toy . . 5
## 1054 -ss SS SING;CONT the brother taste -ss the tea -s . . 4
## 1301 the DET PLUR;CONT the clerk -s drink the tea . . 5showWordOutError("the,a,X.ss,X.ed",c("Agreement","measure","coul"),postverbcoul.df)
corrHiddenOutputInputDeriv(postverbcoul.df,"Probability", "Agreement")## Condition
## ACTCONT CATCONT CATVIOL
## 0 0 0
## CONG GARDAMB GARDSC
## 0 0 0
## HIGHCLOZE INCOH LOWCLOZE
## 0 0 0
## PASSCONT PLUR;CONT PLUR;VIOL
## 0 14 2
## SEMATTR SING;CONT SING;VIOL
## 0 12 4
## STRONG STRONGCONS;EXPECTED STRONGCONS;UNEXPECTED
## 0 0 0
## SUBCONT SUBVIOL TENSECONT
## 0 0 0
## TENSEVIOL TEST-PRED TEST-UNPRED
## 0 0 0
## WEAK WEAKCONS;EXPECTED WEAKCONS;UNEXPECTED
## 0 0 0
## ́ZEROCLOZE
## 0
# make erp figure with different measures
#mapping = aes(x=depth, y=value, colour = Agreement,linetype=Probability)
coulnorawpostverbnotargout.df=removeWordsDrawERP(postverbcoul.df,"value","Probability","Agreement")## Probability Number Agreement sub coul devellearn epoch word
## 1 Probable Singular2 Control s0 0.1vcoul1 <NA> 100000 -ss
## 2 Probable Singular2 Control s0 0.1vcoul1 <NA> 100000 -ss
## 3 Probable Singular2 Control s0 0.1vcoul1 <NA> 100000 -ss
## 4 Probable Singular2 Control s0 0.1vcoul1 <NA> 100000 -ss
## 5 Probable Singular2 Control s0 0.1vcoul1 <NA> 100000 -ss
## 6 Probable Singular2 Control s0 0.1vcoul1 <NA> 100000 -ss
## wordcat contfunc Example Condition tick layer
## 1 SS F the boy take -ss a stick . . SING;CONT 4 word
## 2 SS F the boy take -ss a stick . . SING;CONT 4 word
## 3 SS F the boy take -ss a stick . . SING;CONT 4 compress
## 4 SS F the boy take -ss a stick . . SING;CONT 4 compress
## 5 SS F the boy take -ss a stick . . SING;CONT 4 hidden
## 6 SS F the boy take -ss a stick . . SING;CONT 4 hidden
## measure depth variable value
## 1 outderiv 1 abssum 2.9467200
## 2 inputderiv 1 abssum 1.3212796
## 3 outderiv 2 abssum 22.4376145
## 4 inputderiv 2 abssum 0.2291372
## 5 outderiv 3 abssum 5.0740166
## 6 inputderiv 3 abssum 0.3905577interactionText="NO TEXT"
modellist2 = list()
# get example for table
agreeinputword = subset(coulnorawpostverbnotargout.df,measure == mainMeasure & layer == "word")
exampleTable = agreeinputword[1:4,c("Number","Agreement","Example")]
tabletext = kable(exampleTable,format="latex",row.names=FALSE)
print(exampleTable)## Number Agreement Example
## 159 Singular2 Control the boy take -ss a stick . .
## 393 Singular1 Control the father -s taste wine -s . .
## 614 Singular2 Violation the mother push a kite . .
## 848 Singular1 Control the wife -s take the toy . .coulagreeinput = subset(coulnorawpostverbnotargout.df,measure == mainMeasure & layer %in% c("word",p600layer) )
if (nlevels(adultdata.df$sub) > 4){ # do mixed with more than 9 subj
coulagreeinput$cplural = ifelse(coulagreeinput$Number=="Plural",0.5,-0.5)
coulagreeinput$cviolation = ifelse(coulagreeinput$Agreement=="Violation",0.5,-0.5)
coulagreeinput$cp6layer = ifelse(coulagreeinput$layer==p600layer,0.5,-0.5)
coulagreeinput$cprob = ifelse(coulagreeinput$Probability=="Probable",0.5,-0.5)
#- cviolation:cp6layer:cprob -cviolation:cp6layer - cviolation:cprob - cp6layer:cprob - cp6layer
# this is the main mixed model with centered variables
randomeff = paste(" + (1 + cviolation + cprob |sub)") #Agreement * layer
if (noslopes){ randomeff = " + (1 | sub)"}
formu = as.formula(paste(dependMeasure,"~ cviolation*cp6layer*cprob",randomeff))
omnimodel = lmer(formu, coulagreeinput,control=lmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=50000)))
print(summary(omnimodel))
modellist2 = modelComparison(omnimodel) # create anova table
# this model is for posthocs
formu2 = as.formula(paste(dependMeasure,"~ Probability*layer",randomeff))
omnimodelFactor = lmer(formu2, coulagreeinput)
# print(summary(omnimodelFactor))
model.lsmobj <- lsmeans(omnimodelFactor, ~ Probability | layer)
posthocs = summary(as.glht(pairs(model.lsmobj)))
print(posthocs)
interactionText = printInteraction(posthocs)
}## Linear mixed model fit by REML ['lmerMod']
## Formula:
## abssum ~ cviolation * cp6layer * cprob + (1 + cviolation + cprob |
## sub)
## Data: coulagreeinput
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 50000))
##
## REML criterion at convergence: 5029.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1148 -0.3718 -0.0235 0.1309 6.6473
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.0988485 0.31440
## cviolation 0.0003622 0.01903 -1.00
## cprob 0.2684028 0.51808 -1.00 1.00
## Residual 3.8215845 1.95489
## Number of obs: 1200, groups: sub, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.8428 0.1219 23.320
## cviolation 1.8140 0.1412 12.846
## cp6layer 2.3738 0.1411 16.825
## cprob -1.2622 0.2162 -5.838
## cviolation:cp6layer 2.6469 0.2822 9.381
## cviolation:cprob -0.1361 0.2822 -0.482
## cp6layer:cprob -2.0513 0.2822 -7.270
## cviolation:cp6layer:cprob -0.5235 0.5643 -0.928
##
## Correlation of Fixed Effects:
## (Intr) cviltn cp6lyr cprob cvlt:6 cvltn: cp6ly:
## cviolation -0.035
## cp6layer 0.000 0.000
## cprob -0.845 0.032 0.000
## cvltn:cp6ly 0.000 0.000 0.000 0.000
## cviltn:cprb 0.000 -0.599 0.000 0.000 0.000
## cp6lyr:cprb 0.000 0.000 -0.600 0.000 0.000 0.000
## cvltn:cp6l: 0.000 0.000 0.000 0.000 -0.600 0.000 0.000
## [1] "remove . ~ . - cviolation:cp6layer:cprob "## refitting model(s) with ML (instead of REML)## Data: coulagreeinput
## Models:
## model2: abssum ~ cviolation + cp6layer + cprob + (1 + cviolation + cprob |
## model2: sub) + cviolation:cp6layer + cviolation:cprob + cp6layer:cprob
## model: abssum ~ cviolation * cp6layer * cprob + (1 + cviolation + cprob |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 14 5046.5 5117.8 -2509.2 5018.5
## model 15 5047.6 5124.0 -2508.8 5017.6 0.8653 1 0.3523
## [1] "########## Above comparison for cviolation:cp6layer:cprob "
## [1] "remove . ~ . - cp6layer:cprob "## refitting model(s) with ML (instead of REML)## Data: coulagreeinput
## Models:
## model2: abssum ~ cviolation + cp6layer + cprob + (1 + cviolation + cprob |
## model2: sub) + cviolation:cp6layer + cviolation:cprob
## model: abssum ~ cviolation + cp6layer + cprob + (1 + cviolation + cprob |
## model: sub) + cviolation:cp6layer + cviolation:cprob + cp6layer:cprob
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 13 5096.5 5162.6 -2535.2 5070.5
## model 14 5046.5 5117.8 -2509.2 5018.5 51.973 1 5.627e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cp6layer:cprob *** "
## [1] "remove . ~ . - cviolation:cprob "## refitting model(s) with ML (instead of REML)## Data: coulagreeinput
## Models:
## model2: abssum ~ cviolation + cp6layer + cprob + (1 + cviolation + cprob |
## model2: sub) + cviolation:cp6layer
## model: abssum ~ cviolation + cp6layer + cprob + (1 + cviolation + cprob |
## model: sub) + cviolation:cp6layer + cviolation:cprob
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 12 5094.7 5155.8 -2535.3 5070.7
## model 13 5096.5 5162.6 -2535.2 5070.5 0.2238 1 0.6362
## [1] "########## Above comparison for cviolation:cprob "
## [1] "remove . ~ . - cviolation:cp6layer "## refitting model(s) with ML (instead of REML)## Data: coulagreeinput
## Models:
## model2: abssum ~ cviolation + cp6layer + cprob + (1 + cviolation + cprob |
## model2: sub)
## model: abssum ~ cviolation + cp6layer + cprob + (1 + cviolation + cprob |
## model: sub) + cviolation:cp6layer
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 11 5204.3 5260.3 -2591.2 5182.3
## model 12 5094.7 5155.8 -2535.3 5070.7 111.64 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation:cp6layer *** "
## [1] "remove . ~ . - cprob "## refitting model(s) with ML (instead of REML)## Data: coulagreeinput
## Models:
## model2: abssum ~ cviolation + cp6layer + (1 + cviolation + cprob | sub)
## model: abssum ~ cviolation + cp6layer + cprob + (1 + cviolation + cprob |
## model: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 10 5218.0 5268.9 -2599.0 5198.0
## model 11 5204.3 5260.3 -2591.2 5182.3 15.68 1 7.501e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cprob *** "
## [1] "remove . ~ . - cp6layer "## refitting model(s) with ML (instead of REML)## Data: coulagreeinput
## Models:
## model2: abssum ~ cviolation + (1 + cviolation + cprob | sub)
## model: abssum ~ cviolation + cp6layer + (1 + cviolation + cprob | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 9 5411.5 5457.3 -2696.7 5393.5
## model 10 5218.0 5268.9 -2599.0 5198.0 195.44 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cp6layer *** "
## [1] "remove . ~ . - cviolation "## refitting model(s) with ML (instead of REML)## Data: coulagreeinput
## Models:
## model2: abssum ~ (1 + cviolation + cprob | sub)
## model: abssum ~ cviolation + (1 + cviolation + cprob | sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2 8 5432.2 5472.9 -2708.1 5416.2
## model 9 5411.5 5457.3 -2696.7 5393.5 22.738 1 1.857e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "########## Above comparison for cviolation *** "
## $`layer = word`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Improbable - Probable == 0 1.0470 0.2633 3.977 6.99e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
##
## $`layer = hidden`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Improbable - Probable == 0 3.0983 0.2633 11.77 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)xtabs(~ coul + Agreement, coulagreeinput)## Agreement
## coul Control Violation
## 0.1vcoul1 480 120
## 0.1vcoul2 120 480xtabs(~ Probability + Agreement, coulagreeinput)## Agreement
## Probability Control Violation
## Improbable 120 120
## Probable 480 480# mapping = aes_string(x="depth", y=dependMeasure, colour = "Probability",linetype="Probability")
p = drawERP(coulagreeinput,dependMeasure,"Probability",iv4="Agreement", span=1, timeline.y=4,timesize=2,xmax=7)
# p = p + facet_wrap(~ Agreement)
p
if (!noslopes){ # when doing final version, change noslopes to FALSE
myggsave("img/agreelearn.png",width=figwidth,height=figheight)
}CUT AND PASTE INTO LATEX DOC click in box and use right arrow to see whole text
```latex
This is an example of the stimuli used in the model.
The maximal model for the data had random slopes + (1 + cviolation + cprob |sub)
Mixed Model Results a main effect of cviolation (\(\beta\)=1.8, SE=0.14, \(\chi^2\)(1)=22.74, p<0.001) a main effect of cp6layer (\(\beta\)=2.4, SE=0.14, \(\chi^2\)(1)=195.44, p<0.001) a main effect of cprob (\(\beta\)=-1.3, SE=0.22, \(\chi^2\)(1)=15.68, p<0.001) an interaction of cviolation:cp6layer (\(\beta\)=2.6, SE=0.28, \(\chi^2\)(1)=111.64, p<0.001) no interaction of cviolation:cprob (\(\beta\)=-0.14, SE=0.28, \(\chi^2\)(1)=0.22, p=0.636) an interaction of cp6layer:cprob (\(\beta\)=-2.1, SE=0.28, \(\chi^2\)(1)=51.97, p<0.001) @@no interaction of cviolation:cp6layer:cprob (\(\beta\)=-0.52, SE=0.56, \(\chi^2\)(1)=0.87, p=0.352) Posthoc tests for each layer @@Posthoc: There was a significant difference for layer = word, diff = 1.047, t(0)=3.98, p<0.001 @@Posthoc: There was a significant difference for layer = hidden, diff = 3.0983, t(0)=11.77, p<0.001 ratio layer = hidden 2.95921680993314 rev 0.337927250427654